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Abstract 


The  FORTRAN  program  SCAT  is  used  to  solve  the  magnetic  -field 
integral  equation  (MFIE)  to  determine  the  -fields  scattered  by  a 
perfectly  conducting  sphere.  The  incident  field  is  a plane-wave 
pulse,  and  a steppi ng-i n-ti me  procedure  is  used  to  determirie  the 
surface  current  density  induced  on  the  sphere.  The  program  does 
not  take  advantage  of  the  special  symmetry  of  the  scatterer 
because  it  is  intended  to  serve  as  a verified  starting  point  for 
more  general  programs.  The  output  is  compared  to  that  of  the 
program  PERF,  which  computes  the  same  fields  via  a Fourier 
transform  of  the  monochromatic  fields  obtained  from  the  Mie 
formulas.  The  contributions  of  the  self-patch  and  neighboring 
patches  to  the  singular  integral  are  optionally  computed  by  using 
their  expansions  in  the  linear  size  of  the  patches.  The  self- 
patch term  is  important  for  the  solution  of  other  integral 
equations  that  may  be  of  the  first  kind.  For  the  MFIE,  these 
corrections  are  small  but  not  negligible.  The  program  takes 
advantage  of  the  vector  programming  features  of  the  CYBER  205. 

Key  words:  computer  programs;  el ectromagnet i c pulse  scattering; 

magnetic  field  integral  equation;  Mie  scattering; 
perfectly  conducting  sphere;  stepping-in-time 
procedure;  transient  el ectromagnet i c fields;  vector 
pr ogrammi ng . 
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There  are  many  procedures  related  to  the  metrology  mission 

of  the  National  Bureau  of  Standards  that  require  the  analysis  of 

waves  scattered  by  an  obstacle.  Comparison  of  the  scattered 

waves  with  those  computed  for  known  obstacles  can  serve  to 

determine  parameters  related  to  the  scatterer.  Codes  used  in 

inverse  scattering  problems  also  may  require  the  solution  of  the 

direct  problem.  If  the  scatterer  is  a homogeneous  body,  the 

scattered  fields  can  be  computed  from  tangential  fields  defined 

on  the  surface  of  the  scatterer.  In  turn,  these  fields  can  be 

computed  by  solving  a singular  integral  equation. 

The  el ectromagneti c fields  scattered  by  a perfect  conductor 

* 

can  be  obtained  form  the  surface  current  density  J CID  . This 

s 

tangential  field  can  be  determined  by  solving  the  magnetic  field 
integral  equation  (MFIE) , 

-+  i n 

J (x,t)  = (2/M,-)n  X B (x,t) 


^ n X 4)  dS^R 
2tt 


1 S ■ 1 


(1) 


-♦in 

where  B is  the  magnetic  field  of  the  incident  pulse,  Mq  is  the 
permeability  of  free  space,  n is  the  outside  normal  to  the 


Numbers  in  square  brackets  indicate  the  literature  references 
at  the  end  of  his  report.  There  are  further  references  in  these 
publ i cat i one. 
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scatterer,  t is  the  retarded  time 


T = t - R/c, 


(2) 


and 

R = X - X' , R = IRI  . (3) 


This  is  a singular  integral  equation  because  the  integrand 
diverges  Mhere  R =0,  and  it  is  an  integral  equation  o+  the 
second  kind  because  the  unknown  -field  also  appears  outside  the 
integral,  on  the  le-ft  side  of  the  equation. 

An  example  of  a singular  integral  equation  of  the  first  kind 
is  the  electric  field  integral  equation  (EFIE) , 


^in  -♦  1 " 

0 = nxE  (x,t)  + rr-f 

4n 


X dS' 
S 


(x ' , T ) 


+ 


(x  ' , T ) 

R'^  ® 


Mq  ,T)' 

R aV 


(4) 


1 n 

where  E is  the  electric  field  of  the  incident  pulse,  is  the 

permittivity  of  free  space,  and  is  the  surface  charge  density, 
related  to  by  conservation  of  charge. 

We  have  shown  C2,  3D  how  the  scattered  and  transmitted 
fields  of  a pulse  scattered  off  a homogeneous  dielectric  body  can 
be  determined  from  a single  tangential  vector  field  that  obeys  a 
singular  integral  equation  of  the  first  kind.  We  have  further 
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generalized  C41  this  -formulation  to  the  scattering  by  a 
dispersive  body,  such  as  an  imper-fect  conductor.  In  some  papers 
we  also  derive  equations  -for  the  scattering  of  other  types  of 
fields  by  homogeneous  bodies,  as  well  as  the  correspondi ng 
equations  for  the  scattering  of  monochromat i c waves. 

The  solution  of  the  integral  equations  by  means  of  the 
stepping-in-time  procedure  takes  advantage  of  causality  and  the 
finite  speed  of  propagation  of  the  wave.  The  new  values  of  the 
unknown  field  are  expressed  in  terms  of  an  integral  over  the 
retarded  values  of  the  field,  which  have  been  calculated  in 
previous  steps. 

The  integration  can  be  performed  numerically  by  dividing  the 
surface  of  the  scatterer  into  patches  and  adding  the  products  of 
the  integrand  at  the  center  of  each  patch  times  the  area  of  the 
patch.  This  prescription  does  not  apply  to  the  contribution  of 
the  self -patch  because  this  is  the  patch  where  the  integrand 
becomes  infinite  as  the  field  point  x and  the  source  point  x' 
coincide.  For  an  integral  equation  of  the  second  kind,  the 
contribution  of  the  self— patch  can  be  neglected  in  the 
computation  of  a new  value  of  the  field  in  comparison  to  the 
large  number  of  contr i but i ons  from  the  other  patches.  This  is 
not  so  for  an  integral  equation  of  the  first  kind,  because  the 
self-patch  contribution  to  the  integral  is  the  only  term  that 
depends  on  the  new  value  of  the  field. 

We  have  expanded  the  contribution  of  the  self-patch  to  first 
order  in  the  linear  size  of  the  patch  for  an  arbitrary  surface 
[51,  and  in  this  form  it  may  be  included  in  the  integral.  We 


have  similarly  expanded  the  contribution  o-f  a neighboring  patch, 

because  the  integrand  cannot  be  well  approximated  by  its  value  at 

2 

the  center  o-f  such  a patch  when  1/R  or  1/R  varies  significantly 
over  the  patch  C63. 

Me  have  developed  the  program  SCAT  to  implement  the 
computation  of  the  surface  current  density  and  the  far  fields 
scattered  by  a perfectly  conducting  sphere  by  solving  the  MFIE. 
In  this  program  we  provide  the  option  of  computing  the  correction 
to  the  surface  fields  due  to  the  self— patch  contribution  and  the 
more  accurate  computation  of  the  nei ghboring-patch  contributions. 

We  found  that  terms  that  contain  the  time  derivatives  of  the 
fields,  terms  that  are  negligible  in  most  computations,  cannot  be 
ignored  when  the  fields  themselves  are  small  C73.  The  expansions 

-4 

of  the  surface  current  density  J and  its  derivatives  that  we  use 

s 

in  SCAT  are  explained  in  C73;  they  differ  from  those  we  used  in 

£51  and  £61  in  that  the  dependence  of  the  retarded  time  r on  the 

curvilinear  coordinates  u and  v is  taken  into  account  when  a 
partial  derivative  with  respect  to  u or  v is  taken. 

We  intend  SCAT  to  serve  as  a validated  starting  point  for 

other  programs  that  will  determine  the  fields  scattered  by  more 

general  types  of  scatterers.  For  this  reason  we  do  not  take 
advantage  of  special  symmetry  properties  of  the  sphere.  To 
validate  the  program,  we  compare  the  far  fields  to  those  computed 
by  means  of  another  program,  PERF,  which  uses  the  Mie  formulas 
£71  for  monochromati c fields  and  a Fourier  transform.  This 
program  is  discussed  in  section  2. 

The  general  structure  of  SCAT  is  presented  in  section  3. 
Section  4 is  devoted  to  a discussion  of  the  issues  related  to  the 
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computation  of  the  contri buti ons  of  the  patches  neighboring  the 
self-patch,  and  in  section  5 we  examine  the  contribution  of  the 
self-patch  itself.  In  section  6 we  point  out  where  we  use  vector 
programming  as  appropriate  for  the  CYBER  205.  We  describe  a 
number  of  ancillary  programs  in  section  7,  and  we  show  some 
examples  of  results  in  section  8.  Listings  of  computer  programs 
are  given  in  the  appendix. 

We  define  the  incident  plane-wave  pulse  so  that  it  reaches 
the  sphere  at  t = 0.  We  choose  a double-exponential  smoothed  at 
the  beginning  by  a power  factor  for  the  profile  of  the  plane 
wave.  The  incident  magnetic  field  is 

= O,  <5> 

where  0 is  the  unit  step  function,  a is  the  rise  constant  of  the 
pulse,  ^ the  decay  constant,  and 

f=z+a— ct,  (6) 

where  a is  the  radius  of  the  sphere. 

2.  RC99C§[D  EERE 

The  Mie  formulas  give  the  fields  scattered  by  a homogeneous 
sphere  from  an  incident  plane  monochromatic  wave.  For  a perfect 
conductor  and  a wave  incident  along  the  z— axis,  the  scattered  far 
field  components  reduce  to 
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(8) 


where  k is  the  wavenumber,  r is  the  distance  from  the  center  of 
the  sphere  to  point  where  the  field  is  computed,  © is  the  polar 
or  scattering  angle,  0 is  the  azimuthal  or  polarization  angle, 
P^^^  is  the  associated  Legendre  function,  and  a and  b are 
coefficients  given  by 

a (p)  = - j (/!))  /h  (/>:■  , (9) 
n -*n  n 


b (p)  = - Cpj  (p)  D D' , 

n n n 


(10) 
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is  a spherical 


where  j is  a spherical  Bessel  function,  h^^^ 
n n 

Hankel  function  of  the  first  kind,  and 

p = ka  = 2na/X.  <11) 

The  primes  in  eq.  (11)  indicate  derivatives  with  respect  to  p. 

The  Fourier  transform  of  the  incident  wave  (5)  is 

= expCikCz  + a) 1 (^/c) C - ik)“^  - (a  - ik)“^l,  (12) 

w 

where  k = oj/c.  The  coefficient  of  exp(ikz)  multiplies  the 
amplitude  of  the  far  fields  given  in  eqs.  (7)  and  (8).  The 
steady  state  component  of  the  incident  pulse  does  not  contribute 
to  the  scattered  fields. 

In  the  program  PERF,  we  compute  first  the  required  Bessel 

and  Hankel  functions  by  means  of  the  subroutines  BESRJ  and  BESRY 

l83,  then  the  scattering  coefficients  a and  b in  the  subroutine 

n n 

SCACO,  and  the  monochromati c fields  in  the  subroutine  FLD.  The 
inverse  Fourier  transform  is  carried  out  in  subroutine  INUFT  C9] 
or  VINUFT  (its  vector  version);  we  do  not  use  a fast  Fourier 
transform  because  we  determined  that  it  was  better  to  cover  low 
and  high  frequencies  at  different  intervals  than  to  use  an 
equispaced  grid.  The  frequencies  are  selected  in  subroutine 
OMEGA  C91  and  the  intervals  increase  in  a geometric  progression, 
which  gives  a thorough  coverage  of  low  frequencies  and  a sparse 
but  extended  coverage  of  high  frequencies.  The  input  to  the 
program  requires  the  specification  of  a,  ff,  and  the  radius  of  the 


sphere, 

the 

number  of 

scattering  angles 

and  the 

scattering 

angles. 

In 

addition. 

we  have  to  specify 

the  lower 

and  upper 
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^^requency  limits  at  which  the  fields  are  computed  and  the  number 
of  frequency  and  time  points.  The  output  consists  of  the 

intensity  of  the  scattered  fields  for  both  polarizations  of  the 
scattered  fields;  they  are  saved  in  a file  and  they  are  plotted 
using  the  subroutines  DRAW3  C103  or  DRAW4  till. 

3.  SCAT 

In  this  program  we  compute  the  surface  current  density 
induced  on  the  sphere  by  solving  the  MFIE  (1)  in  the  time  domain 
by  a steppi ng— i n-ti me  procedure.  The  scattered  fields  are  then 
obtained  by  integration,  and  we  compute  the  intensities  in  the 
far  zone  for  comparison  with  the  output  of  PERF.  The  magnetic 
field  in  the  radiation  zone  is 


(X,  t) 


4nr 


r X U dS 


aj  (x',T) 

s 

at^ 


(13) 


The  surface  current  density  is  tangential  to  the  surface  and 
has  only  two  independent  components.  To  decrease  the  amount  of 
memory  required  to  run  the  program,  we  decided  to  store  only  the 
components  along  the  meridians  and  parallels  on  the  sphere.  To 
sum  the  contributions  from  different  patches,  we  first  compute 
the  components  of  the  vectors  in  a Cartesian  coordinate  system. 
Once  the  two  tangential  components  of  the  surface  current  density 
are  determined,  they  are  saved  on  disk  for  future  use. 
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Most  o-f  the  execution  time  needed  -for  the  program  is  spent 
computing  the  surface  current  density  by  solving  the  MFIE.  The 
components  of  this  tangential  vector  are  computed  at  the  center 
of  each  patch  starting  at  t = O and  after  each  time  increment  At. 

If  we  ignore  the  self-patch  contribution  to  the  integral  in 
eq.  (1),  the  contributions  from  the  other  patches  can  be  computed 
because  the  fields  at  their  centers  have  previously  been 
calculated  at  the  retarded  times.  We  have  to  choose  a time 
interval  that  is  smaller  than  the  distance  between  any  pair  of 
patches  (divided  by  the  speed  of  light  c)  to  obtain  a retardation 
that  is  sufficiently  large.  It  is  unlikely  that  the  current 
density  will  have  been  calculated  at  precisely  the  retarded  time, 
and  an  interpol ation  is  needed  to  obtain  the  necessary  values. 
The  i nterpol at i on  is  discussed  in  the  next  subsection. 

After  the  values  of  the  surface  current  density  are  computed 
in  this  way,  we  can  perform  the  corrections  for  the  self— patch, 
for  the  neighboring  patches,  or  for  both;  these  corrections  are 
described  in  detail  in  separate  sections. 

3.2  The  Interpol  at  i_gn  Formulas 

The  interpolation  method  C123  turned  out  to  affect 
significantly  the  outcome  of  the  calculations.  The  program  SCAT 
allows  for  two  alternatives.  In  the  first  one,  the  fields  at  the 
retarded  times  are  computed  from  a polynomial  of  degree  4 that 
passes  through  five  points  chosen  so  that  they  fall  on  both  sides 
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o-f  the  retarded  time.  The  relation  has  the  -form 


y(t)  = a^t"^  + + a^t  + 


(14) 


and  we  rescale  the  independent  variable  so  that  it  takes  integer 
values  from  -2  to  2.  In  terms  of  the  correspond! ng  values  of  the 
dependent  variable,  the  coefficients  are 


^4  “ - iy(l)  + iy<0)  - iy(-D  + 54>'<-2>. 


(15) 


a^  = ^y(2)  - + zy(-D  - ^^(-2), 


12' 


12’ 


( 16) 


a„  = - 


24  J 


(2)  + |y(l)  - zy<0)  + |y(-l)  - irv<-2). 


24' 


(17) 


a^  = - Y^<2)  + |y(l)  - §y<-l)  + 


(18) 


a^.j  = y (O)  . 


(19) 


In  the  second  method,  a polynomial  of  degree  2 is  fitted  to  the 

same  five  points  by  a 1 east-squares  approximation.  The 

coefficients  a.  and  a_  are  zero,  and  the  other  coefficients  are 
4 3 


a„  = iy(2)  - Tr/<1>  " kv  <0)  - + ^y(~2), 


14’ 


14 


(20) 
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(21 ) 


= iy<2) 


(1)  - 


a 


O 


35 


12 


y (-1) 


(22) 


The  polynomial  interpolation  gives  more  accurate  results,  but 
leads  more  easily  to  numerical  instabilities  E7,  13].  The  time 
derivatives  of  the  fields  are  calculated  from 

yMt)  = 4a^t"  + 3a^t^  + 2a^t  + a^.  (23) 

The  values  of  the  fields  and  their  derivatives  are  set  to 
zero  for  times  before  the  time  the  incident  pulse  reaches  the 
center  of  the  correspond! ng  patch  on  the  sphere.  The 
interpolated  values  do  not  necessarily  vanish  at  these  times  and 
the  error  increases  when  they  are  used. 

The  interpolation  formulas  are  implemented  in  a form 
suitable  for  the  execution  of  vector  instructions  in  subroutines 
INTRPO  and  INTRPl,  and  INTRP2  and  INTRP3  are  scalar  versions 
needed  only  for  the  computation  of  the  derivatives  of  the  fields 
at  the  self— patch.  The  choice  of  the  five  points  for  the 
interpolation  also  has  an  effect  on  the  accuracy  of  the  results, 
which  is  improved  if  the  retarded  time  is  close  to  the  center  of 
the  interval.  We  choose  the  points  in  this  way  except  for 
patches  that  are  so  close  to  the  self— patch  that  the  retardation 
is  too  small  for  all  the  needed  values  to  be  available.  The 
values  of  the  components  of  the  surface  current  density  and  of 
their  time  derivatives  at  the  self-patch  itself  are  needed  to 
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carry  out  the  sel-f-  and  neighboring-patch  corrections;  we  use  the 
uncorrected  values  at  the  self-patch  as  a first  approximation 
where  necessary. 

3.3  Subroutine  RUNT 

This  subroutine  is  executed  once  at  the  beginning  of  the 
program  to  perform  the  division  of  the  sphere  into  patches  and  to 
compute  coefficients  that  will  be  needed  during  the  execution  of 
the  program. 

We  do  not  take  advantage  of  the  symmetry  of  the  sphere  to 
simplify  the  integral  equations,  but  we  have  to  take  into  account 
the  character i st i cs  of  the  sphere  to  choose  a good  way  to  divide 
it  into  patches.  We  naturally  use  spherical  coordinates,  the 
polar  angle  0 and  the  azimuthal  angle  which  form  a principal 
orthogonal  curvilinear  coordinate  system  C51. 

We  divide  the  surface  of  the  sphere  into  patches  by 
meridians  and  parallels.  We  first  divide  the  surface  into  n^ 
strips  of  equal  increments  of  0.  Each  strip  is  then  divided  into 
a number  of  patches  that  is  roughly  proportional  to  the  length  of 
the  strip,  that  is,  to  sin0.  We  leave  two  undivided  caps  at  the 
poles,  where  spherical  coordinates  become  undefined;  we  use 
Cartesian  coordinates  for  these  caps,  which  requires  special 
formulas  for  the  computations.  We  divide  the  strips  next  to  the 
caps  into  a given  number  n^  of  patches,  usually  5.  Then  the 
other  strips  are  divided  into  approx i matel y 

n = (n_  - n^)sin0  + n_  (24) 

2 •:>  3 
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patches  o+  equal  size  by  lines  a-f  constant 

This  more  or  less  asymmetrical  division  of  the  surface  into 

patches  leads  to  computed  nonzero  values  for  some  intermediate 

variables  that  should  vanish  by  symmetry.  Although  this 

discrepancy  does  not  seem  to  make  much  difference  in  the  final 

results,  we  get  intermediate  values  that  are  closer  to  zero  by 

restricting  the  number  of  patches  in  each  strip  to  a power  of  2. 

For  strips  near  the  poles,  n is  not  large  and  the  values  of 

the  increment  in  = 2n/n,  may  not  be  small  compared  to  1. 

The  side  of  the  patch,  asinOA^,  is  small  compared  to  a. 

To  obtain  patches  that  are  approx i matel y square  in  form,  we 

would  have  to  choose  n = 2n  . A choice  that  works  essentially 

^ 1 

just  as  well  is  n^  = n^.  The  consequent  reduction  in  memory 

requirements  and  computing  time  is  a made  possible  by  the 
periodicity  of  the  fields  in  0,  which  greatly  improves  the 
accuracy  of  the  integration  over  that  variable  C73. 

We  compute  and  store  the  components  of  the  radius  vector,  of 
the  tangential  vectors,  the  surface  elements,  and  the 

coefficients  required  for  self—  and  nei ghbori ng— patch 

corrections.  These  computations  involve  the  evaluation  of 
tr i gonometr i c functions  and  logarithms,  the  performance  of 
numerical  i ntegrati ons,  and  other  time-consuming  operations,  but 
they  are  done  only  once  in  each  run. 

3.4  Subroutine  INFL 

i n 

The  incident  field  B has  to  be  computed  at  the  center  of 
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this  is  done  in  subroutine  INFL. 


each  patch  at  each  time  step; 

For  the  plane-Mave  pulse,  eq.  (5)  allows  us  to  calculate  the 
magnetic  -field  -for  arbitrary  t.  I-f  the  pulse  is  spatially 

limited,  the  incident  -fields  have  to  be  computed  -from  their 
initial  values  by  solving  Maxwell's  equations.  This  can  be  done 
by  numerical  integration  over  the  information-col lecting  sphere 
C 14]  , 

3.5  Subroutine  FARFL 

The  far  field  is  calculated  in  subroutine  FARFL  by  doing  the 
integration  in  eq.  (13)  for  any  number  of  scattering  directions 
given  by  the  corresponding  values  of  Q and  ^ for  the  direction  of 
the  scattered  pulse.  The  intensities  of  the  fields  are  plotted 
and  they  are  also  saved  on  disk  for  comparison  with  the  results 
of  other  cal cul ati ons. 

The  time  scale  is  referred  back  to  the  center  of  the  sphere, 
as  is  appropriate  to  compare  the  fields  with  the  output  from 
PERF.  The  time  the  scattered  pulse  starts  is  a function  of  the 
scattering  angle, 

t.^  = (a/c)Cl  - 2sin<e/2)l,  (25) 

and  the  duration  of  the  pulse  is  comparable  to  the  length  of  the 
pulse  and  to  the  size  of  the  sphere. 
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4. 


Sel-f -Patch  Corrections 


The  surface  integral  in  eq.  (9)  is  singular  because  the 

integrand  diverges  when  R = O,  which  happens  on  the  self-patch. 

The  magnitude  of  the  first  term  in  the  integrand  is  proportional 

to  R ^ and  the  integral  of  that  term  is  well  defined.  The 

-2 

magnitude  of  the  second  term  is  proportional  to  R , and  the 
integral  is  defined  in  a limiting  sense  analogous  to  the 
principal  value  of  a one-dimensional  integral.  In  an  expansion 
in  the  size  of  the  patch  C53,  the  leading  term  in  the  integral  is 
not  absolutely  integrable,  but  it  vanishes  by  symmetry  when  the 
field  is  computed  at  the  center  of  the  patch.  The  next  term  is 
linear  in  the  size  of  the  patch. 

The  self-patch  contribution  to  the  integral  is  usually 
neglected,  although  there  is  no  more  reason  to  neglect  the 
contribution  of  this  patch  than  that  of  any  other  patch. 
Actually,  the  self— patch  contribution  is  linear  in  the  size  of 
the  patch,  while  those  of  most  other  terms  are  quadratic.  Also 
the  seif— patch  contribution  is  of  fundamental  importance  in  an 
integral  equation  of  the  first  kind,  where  the  unknown  field  does 
not  appear  outside  the  integral.  We  thus  include  in  the  program 
the  option  to  compute  the  self— patch  correction. 

The  expansion  of  the  integral  of  the  first  term  is  of  order 
2 in  the  size  of  the  patch,  but  we  cannot  neglect  this  term  when 
is  comparable  to  L73.  In  this  case,  because  cAt  is  less 
than  R in  the  denominators,  the  integral  of  the  first  term  is 
larger  than  the  integral  of  the  second  term.  This  happens,  for 
instance,  when  the  pulse  first  reaches  a patch.  To  make  sure 
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we  keep  this 


that  we  do  not  neglect  an  important  contribution, 
term  in  all  self-patch  corrections.  The  importance  of  this  term 
complicates  the  implementation  of  a similar  solution  for  an 
integral  equation  of  the  first  kind. 

We  include  terms  with  spatial  derivatives  of  the  surface 
current  density,  which  are  often  ignored  or  assumed  to  be 
negligible  by  requiring  that  the  fields  vary  slowly  over  the 
scatterer  H153.  We  also  have  to  take  into  account  the  spatial 
variation  of  the  tangential  unit  vectors. 

The  resulting  expression  for  the  self  patch  correction  is 
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where  the  coefficients  are  given  in  C51  or  C71,  0 and  ^ are  thi 
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tangential  unit  vectors  along  the  coordinate  curves,  and  A9  is 
the  angular  width  of  a strip.  To  compute  the  corrections,  we 
need  to  know  the  fields  and  their  time  derivatives  at  the  self- 
patch, where  the  time  retardation  vanishes.  We  determine  their 
values  from  the  uncorrected  results.  If  the  integral  of  the  term 
proportional  to  SJ^/dt  is  neglected,  the  correction  is 

proportional  to  and  it  can  be  incorporated  in  the  coefficient 
of  the  term  on  the  left  side  of  the  equation. 

Derivatives  with  respect  to  <6  are  calculated  by 

approximating  them  by  a quotient  of  differences  of  the  respective 
quantities  for  the  patches  on  both  sides.  Derivatives  with 
respect  to  0 are  calculated  in  a similar  manner.  If  the  values 
of  the  functions  at  the  same  ^ on  the  strips  above  and  below  are 
not  available,  they  have  to  be  obtained  by  interpolation.  There 
IS  an  additional  problem  when  one  of  the  neighboring  strips  is  a 
polar  cap  because  the  0-  and  fli -components  of  a tangential  vector 
are  not  well  defined  there?  we  define  the  0-direction  tangent  to 
the  meridian  that  passes  through  the  center  of  the  self-patch, 
and  the  correspondi ng  tfk-direction  from  = n x 0.  At  the  poles, 
the  derivatives  with  respect  to  x and  y are  obtained  from  the 
values  on  both  sides  on  the  neighboring  strip. 

Since  there  is  only  one  self-patch  contribution  to  the 
integral,  the  effect  is  hardly  noticeable. 

5.  Neighboring-Patch  Cgrrectigns 

The  integral  over  patches  other  than  the  self— patch  is  first 
computed  by  multiplying  the  value  of  the  integrand  at  the  center 
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o-f  the  patch  by  the  sur-face  o-f  the  patch.  This  is  a reasonable 

approach  it  the  integrand  varies  slowly  over  the  patch.  This  is 

2 

not  so  for  factors  like  1/R  or  1/R  where  R is  small,  that  is, 
for  patches  near  the  self-patch.  Neighboring  patches  are  those 
that  are  closer  than  a small  multiple  (set  in  subroutine  RUNT)  of 
the  patch  size  to  the  center  of  the  self-patch.  We  expand  the 
contributi ons  from  the  neighboring  patches  in  the  parameters  that 
are  of  the  order  of  magnitude  of  the  size  of  the  patch,  which 
here  includes  the  displacements  between  the  centers  of  the 
patches  C63. 

As  we  explained  in  the  previous  section,  the  contribution  of 
the  time-derivative  term  is  not  negligible  when  the  current 
density  at  the  previous  time  is  zero  or  small,  and  we  also  have 
to  expand  this  term  L71. 

We  have  to  compute  n x I,  where 
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tangential.  Also  the  change  in  n between  neighboring  patches  is 


not  small  when  these  patches  are  on  a strip  near  a pole  of  the 
sphere.  Consequently,  we  expand  I itself  instead  of  n x I . We 
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where  6'  and  are  the  tangential  unit  vectors  at  x",  the  center 
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of  S^,  jo'  is  the  radial  unit  vector  in  the  xy-plane,  and  the 
current  density  and  its  derivatives  are  evaluated  at  at  the 

correspond! nq  retards  time  The  derivatives  with  respect  to  0 

or  ^ take  into  account  the  dependence  of  the  retarded  time  on  0 
and  The  coefficients  are  given  in  C63  and  C71. 

Some  of  the  coefficients  become  undefined  when  one  of  the 
sides  of  the  neighboring  patch  has  the  same  coordinate  as  the 
center  of  the  sef-patch  C73.  We  then  have  to  evaluate  these 
coefficients  from  alternative  expressions.  Also  the  azimuthal 
angles  of  the  first  and  the  last  patches  on  a strip  differ  by  an 
amount  close  to  2nj  we  have  to  subtract  2n  from  such  a difference 
to  get  the  proper  value  of  for  the  calculations. 

If  the  neighboring  patch  is  the  spherical  cap  at  0 = O, 
n,  we  use  Cartesian  coordinates  and  obtain 

I it  (0,0,3)  - iC'  <1,0,3)  - jCM0,l,3) 

+ <0,0,3)  - C'  <2,0,3)  - CM0,2,3)jy  x 

+ []o^jCMl,0,3)  - iCM2,0,3)  - jC'<l,l,3)j  x 3J^/dx' 

+ p^CM0,l,3)  - iCMl,l,3)  - jCM0,2,3)j  x dJ^/dy' 
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+ ^ -(fl^CM0,0,2)  - iCMl,0,2)  - jCM0,l,2) 

^ “ CM2, 0,2)  -CM0,2,2)j^  x 3J^/3t' 

+ - f2.CMl,0,2)  - iCM2,0,2)  - jCMl, 1,2)1  x 3J  /3t' 

cLO  J 3x'  s 

+ - f2^CMO,l,2)  - iCM  1,1,2)  - j’CM0,2,2)l  x /3t',  (29) 

clO  ’ -*  J 3 y s 

where 

-f  ■*’' 

= asine  , jsin*^),  (30) 

and  the  coefficients  are  defined  in  C7J.  The  computation  of 
these  coefficient  involves  a numerical  integration  C7D,  for  which 
we  use  the  subprogram  QIDB  in  CMLIB  C161.  The  integrand  becomes 
undefined  for  certain  values  of  the  variable  of  integration,  and 
an  alternative  expression  has  to  be  used  where  this  happens. 

The  corrected  nei ghboring-patch  contr i but i ons  replace  those 
obtained  in  the  simple  way.  We  have  compared  the  two  computed 
values  of  the  nei ghbor i ng-patch  contributions  and  we  have  found 
that  sometimes  they  differ  little  while  at  other  times  the 
difference  can  be  large.  Since  a patch  can  have  about  ten 
neighboring  patches,  the  overall  effect  of  this  correction  is 
more  significant  than  that  of  the  self-patch. 
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6„  Vector  PrQacamming  Considerations 

The  largest  savings  in  execution  time  can  be  achieved  by 
vectorizing  the  computations  in  the  main  loop. 

The  values  of  the  components  of  that  may  be  needed  for 
the  integration  are  kept  in  memory  in  two-dimensional  arrays, 
where  one  index  represents  the  patch  number  and  the  other  the 
time.  After  these  fields  are  no  longer  needed,  they  are  saved  on 
a disk  file. 

An  important  observation  that  facilitates  the  vectori zat i on 
of  the  i nterpol at i on  in  the  time  variable  is  that  the  location  of 
the  retarded  fields  in  the  two-dimensional  array  for  a fixed 
field  point  is  the  same  for  the  different  components  of  the 
fields  and  is  shifted  by  one  column  as  the  time  is  incremented 
for  the  interpol ation.  Thus  the  array  IND  is  computed  once  and 
used  in  a QSVGATHR  function  repeatedly.  This  function  is  used  to 
collect  the  values  of  the  components  of  at  the  different 
retarded  time  in  a single  array.  We  originally  computed  the 
values  stored  in  the  array  IND  and  the  related  array  FT  once  for 
each  field  point  and  saved  them,  to  be  used  at  each  time.  This 
procedure  requires  square  arrays  for  IND  and  FT,  increasing  the 
memory  requi rements.  Although  the  CYBER  205  has  virtual  memory, 
paging  generates  much  activity  and  slows  significantly  the 
execution  of  the  program.  The  arrays  that  contain  the  values  of 
the  surface  current  density  components  are  accessed  in  an 
irregular  manner  and  they  have  to  be  in  memory  at  the  time  of 
execution.  The  arrays  IND  and  FT  may  be  paged  in  gradually  as 


the  field  point  changes,  but  we  have  chosen  to  recompute  their 


values  -for  each  time  point  with  only  a small  increase  in  the 
running  time.  We  save  square  bit  arrays  BT3  and  BT4  to  control 
where  the  interpolation  has  to  be  shi-fted  from  the  center  of  the 
i nterval . 

The  vector i zati on  of  this  main  loop  produces  a reduction  of 
the  execution  time  by  a factor  of  7 if  the  sphere  is  divided  into 
70  patches,  and  20  for  550  patches.  These  factors  can  probably 
be  reduced  somewhat  by  improving  the  scalar  version  of  the 
program,  where  we  use  the  three  Cartesian  components  of  the 
fields  instead  of  the  two  tangential  components  we  use  in  the 
vector  program.  The  vectori zati on  has  to  be  done  explicitly  to 
be  able  to  take  advantage  of  the  Q8VGATHR  function. 

The  subroutine  RUNT  is  called  only  once  at  the  beginning  of 
the  program,  and  we  made  no  great  effort  to  vectorize  this  part 
of  the  code.  The  subroutine  FARFL  can  be  called  repeatedly  to 
find  the  scattered  fields  in  different  directions,  and  its 
structure  is  not  very  different  from  the  main  loop  since  we  are 
also  doing  surface  integrals  over  retarded  fields. 

The  computation  of  the  self— patch  correction  takes  only  an 
additional  one  percent  for  execution.  The  additional  execution 
time  for  the  nei ghbori ng-patch  corrections  depends  on  the  number 
of  neighboring  patches,  which  in  turn  depends  on  the  criteria 
used  to  define  them  and  the  patch  distribution  on  the  sphere.  It 
varies  between  ten  and  40  percent  of  the  basic  execution  time  for 
the  loop.  It  may  be  possible  to  decrease  this  time  by 
vectorizing  the  computation  of  the  nei ghbor i ng-patch 
contributions. 
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7. 


Anci.1  jLary  PCQQCSfDS 


We  have  developed  two  programs  that  we  use  to  examine  the 
results  of  the  main  programs  PERF  and  SCAT.  They  are  PERFPLT, 
and  CURCHK.  These  programs  require  very  little  execution  time, 
but  they  run  on  the  CYBER  205  because  that  is  where  the  output 
-files  are  stored. 

We  have  written  a program  SCATP  that  runs  on  the  front-end 
computer  to  determine  the  required  dimensions  of  the  different 
arrays  in  SCAT  and  the  number  of  large  pages  needed  to  accomodate 
these  arrays. 

We  have  also  written  vector  versions  of  subroutines  NUFT  and 
INUFT,  which  we  use  to  compute  the  direct  and  inverse  Fourier 
transform  of  functions  that  are  not  given  at  equispaced  points. 


7.1  Program  PERFPLT 

Program  PERFPLT  uses  the  output  files  PERFSV  of  PERF  and 
SCFLxxx  of  SCAT,  where  xxx  stands  for  the  qualifier  that 
identifies  the  run,  and  allows  us  to  plot  the  far— field 
intensities  on  a single  plot  for  each  outgoing  direction  of 
propagation.  We  thus  can  compare  the  plotted  output  of  PERF  with 
up  to  4 plots  produced  by  SCAT.  We  can  choose  printer  plots  by 
using  DRAW3  ClOl  or  plots  that  can  be  produced  on  a Versatec 
pr 1 nter /pi otter  by  using  DRAW4  Cl  ID  adapted  to  FORTRAN  77. 
Actually,  the  information  for  the  plots  is  added  to  a file  stored 
on  the  front-end  computer.  Then,  when  several  plots  have 
accumulated,  this  information  is  transfered  to  tape  which  is 
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taken  to  our  local  PE  3230  minicomputer  -for  plotting. 

Me  can  select  the  initial  and  -Final  times  on  the  plots  to 
show  details  o-F  particular  regions  of  the  pulses. 

7.2  Prograrn  CURCHK 

Program  CURCHK  selects  the  values  of  the  components  of 
saved  in  a file  SCSVxxx  for  up  to  10  patches  for  the  whole  time 
range  or  a given  interval.  These  values  can  be  printed  or  they 
can  be  plotted  on  the  printer  or  on  a plotter  and  they  provide  a 
useful  check  on  the  solution  of  the  integral  equation. 

7.3  Program  SCATP 

We  run  this  program  on  the  CYBER  180/855,  the  front  end,  and 
then  use  the  editor  to  pass  the  information  to  the  parameter 
statements  in  SCAT  and  CURCHK.  The  class  and  priority  of  a job 
on  the  CYBER  205  depends  on  the  memory  and  time  requi rements,  and 
it  is  better  to  keep  them  as  small  as  possible. 

Program  SCATP  does  the  part  of  the  calculations  that  will  be 
made  in  SCAT  to  find  the  number  of  patches  for  the  given  n^,  n^, 

and  n^5  estimates  the  number  of  time  steps  that  have  to  be  kept 
in  memory,  and  the  number  and  distribution  of  neighboring  patches 
if  this  correction  will  be  made. 

A discrepancy  arose  in  these  calculations  due  to  the 
differences  between  the  CYBER  170/855  and  the  CYBER  205.  The 
former  returned  the  result  cos(n/6)  = .4999...  while  the  latter 
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returned  cos(n/6)  = .5j  rounding  led  to  a different  number  of 
patches  and  an  insufficient  size  for  the  arrays;  this  is  why  the 
size  is  incremented  by  2 in  SCATP. 


7.4  Subroutine  VINUFT 


This  subroutine  is  used  by  program  PERF  to  compute  the 

inverse  Fourier  transform  of  the  scattered  fields.  It  is  related 

to  subroutines  NUFT  and  INUFT  [91. 

The  subroutines  NUFT  can  be  used  to  compute  the  Fourier 

transform  of  a function  f (t) , a 5 t 5 b when  the  values  of  f (t) 

are  known  at  a set  of  values  t^  that  are  not  equi spaced.  The 

function  is  approximated  by  straight  line  segments  joining 

adjacent  points,  the  Fourier  transform  of  this  approx i mat i on  is 

determined  anal yti cal  1 y , and  the  result  is  evaluated  numerically 

at  a given  set  of  values  w. . 

1 

Subroutine  NUFT  can  be  used  twice  to  obtain  the  inverse 
Fourier  transform  in  the  analogous  approximation,  but  we  can  do 
this  more  efficiently  by  using  subroutine  INUFT  if  we  know  that 
f(t)  IS  real,  that  is,  when 

f (-to)  = f *(w:>  . (31 ) 


Subroutines 
NUFT  and  INUFT, 
versions  do  not 
function  used  in 


9NUFT  and  VINUFT  are  the  vectorized  versions  of 
and  they  are  listed  in  the  appendix.  The  vector 
have  the  small  angle  approximation  for  the  sine 
the  scalar  versions. 
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9"  EHQQCSfl!  iQByfes  and  Out  guts 


The  inputs  required  to  run  PERF  and  SCAT  are  contained  in 
the  files  PERFIN  and  SClNxx  (filenames  on  the  front  end  can  have 
7 characters  only,  and  we  add  another  digit  to  indicate  the  type 
of  correction  selected).  The  procedure  used  to  submit  a run  to 
the  computer  causes  a file  SCINxxx  to  be  saved,  so  that  the  input 
file  SCINxk  may  be  modified. 

The  information  that  is  required  for  a given  run  relates 
both  to  the  physical  problem  and  to  the  numerical  computation. 

The  constants  a and  ff  determine  the  shape  of  the  pulse,  and 
we  also  need  the  radius  a of  the  sphere.  Actually,  the  value  of 
a simply  fixes  the  scale  of  the  problem,  and  we  could  set  a = I 
everywhere  in  the  program.  The  outputs  for  a = 1,  a,  B,  and  t 
are  the  same  as  those  for  a,  a/a,  fi/a,  and  at. 

We  have  to  specify  the  angles  0 and  ^ of  the  direction  of 
propagation  of  the  scattered  fields,  and  we  may  give  times  t^  and 
t„  to  obtain  plots  of  some  part  of  the  scattered  fields.  Several 
sets  of  directions  and  times  may  be  specified,  and  all  the 
cai cui at i ons  use  the  same  previously  determined  surface  currents. 
If  tj^  or  t^  is  not  given,  we  make  an  estimate  of  the  extreme 
values  based  on  the  time  of  arrival  of  the  pulse  and  the  time 
covered  by  the  calculation  of  J^.  Further  sets  of  fields  may  be 
computed  in  separate  runs  without  recomputing  the  surface  current 
density,  which  is  saved  in  SCSVxxx. 

The  parameters  needed  for  the  computation  of  the  fields  in 
SCAT  are  the  time  increment  At  and  the  patch  distribution  on  the 
sphere.  The  distance  light  travels  in  the  time  At  has  to  be 
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smaller  than  the  separation  betMeen  the  centers  of  any  two 
patches;  the  program  SCAT  makes  this  At  equal  to  O.S  times  the 
smallest  separation  between  the  centers  if  none  is  given  or  if 
the  given  At  is  larger  than  this  quantity.  We  can  also  give  the 
number  of  time  increments  in  the  calculation,  which  determines 
trie  time  interval  for  which  the  surface  current  density  is 
determined;  if  this  number  is  not  given,  the  program  selects  a 
number  that  corresponds  to  propagation  over  a distance  6a.  The 
division  of  the  sphere  into  patches  is  controlled  by  the  number 
of  strips,  n^,  the  approximate  number  of  patches  on  the 

equatorial  strip,  n^,  and  the  number  of  patches  on  the  strip  next 
to  a polar  cap,  n^,  as  discussed  in  section  3.  To  do  the 

calculation  is  PERF  we  also  have  to  give  the  first  frequency 
increment,  the  highest  frequency,  and  the  number  of  frequency 
values  where  the  monochromat i c fields  are  to  be  calculated. 

The  agreement  between  the  curves  computed  via  the  Fourier 

transform  and  by  solving  the  MFIE  is  shown  in  figure  1.  The 

output  of  SCAT  in  this  figure  corresponds  to  the  largest  run  we 

made.  We  needed  2S  large  pages,  which  is  essentially  the  size  of 

the  computer  memory,  the  main  loop  executed  in  10754  s with  no 

self-  or  nei ghbori ng— patch  corrections,  and  the  total  CPU— time 

for  the  job  was  10812  s.  The  sphere  was  divided  into  58  strips, 

and  this  was  also  the  maximum  specified  for  the  number  of  patches 

in  a strip  (in  powers  of  2,  the  actual  number  goes  from  8 near 

the  poles  to  64  at  the  equator),  for  a total  of  2898  patches.  We 

-10 

selected  At  = 0.3  x 10  s for  the  time  interval  and  computed 
the  current  density  for  900  time  steps.  There  is  some  evidence 
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SCATTERED  INTENSITY,  HORIZONTAL  POLARIZATION  SCATTERED  INTENSITY,  HORIZONTAL  POLARIZATION 


Intensity  (in 
i a -Function  of 


Figure  1. 
fields  a 

and  (c)  at  90  degrees  with  both  pol arizati ons, 
backward  direction  for  an  incident  pulse  with  a = 


f ar 
(b) 
the 
* 

m . The  solid  curves  are  computed  via  a Fourier  transform 

from  liie  scattering  formulas  and  the  dashed  curves  with  the 
markers  were  computed  via  the  MFIE.  There  is  good  agreement. 


arbitrary  units;  of  the  scattered 
ct  in  (a)  the  forward  direction, 

and_  j ^ 

= '^  a\  and 
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especially  in  backward 


of  oscillations  at  the  later  times, 
scat teri ng . 

Many  factors  affect  the  accuracy  of  the  computation  of  the 
fields  from  the  integral  equation,  and  these  factors  interact  in 
complex  ways.  The  two  effects  of  the  i nterpol at i on  formula  used 
for  the  time  variable  are  illustrated  in  figure  2 for  runs  in 
which  the  sphere  was  divided  into  19  strips  with  8 to  32  patches 
in  a strip  (370  patches  total)  or  5 to  20  patches  in  a strip  (258 
patches  total),  and  with  At  = 2.75  x 10  s,  which  is  a 
relatively  large  interval.  Near  the  first  peak,  the  curves  that 
correspond  to  the  fourth— degree  polynomial  i nterpol at i on  are  much 
closer  to  the  correct  solution,  while  near  the  end  of  the  graph 
these  curves  show  the  onset  of  oscillations.  These  oscillations 
tend  to  become  larger  as  time  progresses,  limiting  the  usefulness 
of  this  interpolation  scheme.  We  also  see  in  figure  2 that  the 
particular  way  in  which  the  strips  are  divided  into  patches, 
whether  we  use  powers  of  2 or  not,  has  little  effect  on  the  final 
r esul ts. 

In  figure  3 we  show  a details  of  the  curves  about  the  first 
peak  to  illustrate  the  effect  of  the  self—  and  nei ghbor i ng-patch 
corrections.  The  sphere  was  divided  into  1986  patches,  and  the 
time  increment  was  0.5  x 10  s.  We  see  that  the  effect  of 
including  the  contribution  of  the  self-patch  is  slight,  while 
that  of  correcting  the  contributions  of  the  neighboring  patches 
is  more  significant. 

The  periodicity  of  the  surface  fields  as  functions  of  the 
azimuthal  angle  ^ has  the  effect  of  reducing  the  number  of 
patches  that  are  required  in  a strip  L71.  The  trapezoidal  rule 
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Figurs  2.  Illustration  of  the  effects  of  the  subdivision  of  the 
strips  and  the  interpolation  scheme.  The  curves  marked  with  O 
and  2 correspond  to  patch  numbers  that  are  powers  of  2 and 
curves  O and  1 are  computed  wit!i  1 east— squares  interpolation. 
The  curves  were  cut  off  at  ct  = 2.8  m. 


gives  a much  more  accurate  integral  for  periodic  than  for 
aroitrary  functions  when  the  integration  is  over  one  period. 

if  the  sphere  is  divided  into  n strips,  it  is  better  to 
divide  the  equatorial  strip  into  n^  patches  instead  of  the  2n^ 
required  to  make  the  patches  more  or  less  square.  The  smallness 
of  the  effect  of  this  or  even  a greater  reduction  is  shown  in 
figure  4.  The  increase  in  the  size  of  the  patch  has  an  adverse 
effect  in  the  computation  of  the  self-  and  nei ghbor i ng-pat ch 
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LO 


Figure  3.  Effect  of  seif-patch  and  neighboring-patch 
corrections.  We  show  a detail  of  the  curves  near  the  first 
maximum.  Curves  2 and  3,  computed  with  nei ghbor i ng-patch 
corrections,  agree  better  with  the  solid  curve  than  0 and  1. 
The  inclusion  of  self— patch  corrections  in  the  computation  of  1 
and  3 does  not  lead  to  a noticeable  improvement  over  O and  2, 
respect i vei y . 


corrections,  which  is  also  illustrated  in  figure  4. 

We  show  the  components  of  the  current  density  on  the  top  cap 
in  figure  5a;  the  x— component  remains  essentially  equal  to  zero, 
as  expected  by  symmetry.  The  computed  surface  current  density 
remains  exactly  equal  to  zero  until  the  incident  wave  reaches  the 
patch;  after  that  time,  the  current  density  is  determined  by  the 
integral  equation.  On  the  other  hand,  the  actual  fields  cannot 
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Figure  4.  1 1 i ustrat i on  o-f  the  small  effect  of  a variation  in  the 

number  of  patches  per  strip.  The  sphere  is  divided  into  29 
strips,  and  the  number  of  patches  in  the  equatorial  strip  is  15 
for  curve  O,  29  for  curve  1,  and  45  for  curve  2.  Curve  3 
includes  corrections  with  the  same  patch  distribution  used 
for  curve  O,  showing  a deter i orat i on  of  the  results. 


reach  a patch  until  they  go  around  the  outside  of  the  perfectly 
conducting  sphere.  This  means  that  the  solution  of  the  integral 
equation  at  the  top  cap  of  the  sphere  has  to  remain  zero  between 
ct  = 2a  and  ct  = (1  + n)a,  an  effect  that  is  illustrated  in 
figures  5b  to  5f . The  di screpanci es  shown  for  the  different 
calculations  are  an  indication  of  the  accuracy  that  was  attained 
in  the  solution  of  the  integral  equation. 

The  fields  computed  for  the  program-sel ected  value  of  At  are 
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Figure  5.  Components  of  the  current  density  (in  arbitrary  units) 
at  the  top  pole  of  the  sphere;  the  x-component  remains  zero. 
We  show  the  complete  graph  (a)  in  the  run  with  2898  patches 
with  1 east-squar es  i nterpol at i on , and  a detail  of  the  region 
between  ct  = 2 m and  ct  = 1 + n/2  m for  (b)  the  same  run,  (c) 
one  with  1586  patches,  (d)  same  with  corrections,  (e)  one  with 
370  patches,  and  (f)  370  patches  with  polynomial  i nter pol at i on . 
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Figure  6.  Fields  computed  with  prog^0m— sel ected  (minimum)  values 
of  At  (approx  i mate!  y 3.5  x 10  s)  and  number  o-f  time 

increments  with  370  patches.  The  computation  without 

corrections  shown  in  curve  O is  compared  to  the  same  with 
polynomial  interpolation  corrections  in  curve  1 and  with  a 

time  increment  o-f  1.  x 10  s in  curve  2. 


compared  with  those  computed  with  a smaller  time  increment  in 
-figure  6.  The  curves  show  that  it  is  important  to  use  the  more 
accurate  polynomial  interpolation  Tor  a small  number  o-f  patches. 

The  scattered  -fields  -for  a sharper  incident  pulse,  which  has 
constants  a = 6m^,  ^ = 9m^,  are  shown  in  -figure  7.  The  fit 
at  the  peaks  is  not  very  good  and  there  are  oscillations  at  later 
times.  Smaller  patches  are  required  for  a sharper  pulse. 
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Figure  7. 

Intensity 
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polarizations,  and 
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pulse  with 

a = 6 

X 

m 

and  ^ = 

correspond  to  runs 

wi  th 

1586 

patches. 

oscillations  after  the  initial  peak. 


far  fields  in  (a)  the 
90  degrees  with  both 
kwar^  direction  for  an 
9 m . The  graphs,  which 
show  an  increase  in  the 
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9. 


Concl usi ans 


The  agreement  between  the  results  of  programs  PERF  and  SCAT 
shows  that  both  programs  are  essentially  correct.  Thus,  program 
SCAT  can  be  used  as  the  basis  for  more  general  programs  that 
compute  scattering  of  pulses  by  arbitrary  perfectly  conducting 
surfaces.  The  parts  of  the  program  that  would  have  to  be 
modified  are  mainly  those  related  to  the  division  of  the  surface 
into  patches  and  the  computation  of  spatial  derivatives.  The 
scattering  of  spatially  localized  pulses  can  also  be  included  in 
the  program  by  making  the  appropriate  changes  in  the  computation 
of  the  incident  fields. 

The  importance  of  the  time-derivative  term  in  the 
computation  of  the  contribution  of  the  self -patch  causes 
difficulties  in  the  use  of  a similar  procedure  to  solve  integral 
equations  of  the  first  kind,  such  as  the  EFIE  or  the  integral 
equation  that  describes  the  scattering  by  a dielectric. 

The  correction  for  the  nei ghbori ng— patch  contributions  has  a 
greater  effect  than  inclusion  of  the  self-patch,  probably  because 
there  are  a number  of  neighboring  patches  for  each  patch.  Gn  the 
other  hand,  the  uncorrected  integrals  also  gave  satisfactory 
solutions  with  a sufficiently  large  number  of  patches. 

The  vector i zati on  of  the  program  reduces  the  running  time  by 
a factor  that  depends  on  the  number  of  patches  on  the  sphere. 
This  factor  is  approx i matel y equal  to  20  for  550  patches.  We 
also  found  that  straight  polynomial  interpolation  gave  better 
accuracy  than  a least-squares  polynomial  interpolation  at  the 
cost  of  increased  instability. 
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AgEendiM 


In  this  appendix  we  provide  listings  -for  the  programs  SCAT, 
SCATP,  PERF,  PERFPLT,  and  CURCHK  together  with  most  of  the 
subprograms  that  they  require.  We  also  provide  listings  for  the 
Fourier  transform  subroutines  VNUFT  and  VINUFT. 

The  plotting  subroutines  DRAW3  and  DRAW4  are  documented 
elsewhere,  as  are  the  integration  subroutine  QIDB,  the  interval 
generating  subroutine  OMEGA,  and  the  subroutines  to  compute 
Bessel  functions  BESRJ  and  BESRY. 
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PROGRAM 


SCAT 


r 

C VECTORIZED  VERSION  OF  PROGRAM  THAT  COMPUTES  SCATTERED  FIELDS. 

C SELF-PATCH  AND  NEIGHBORING-PATCH  CONTRIBUTIONS  CAN  BE  INCLUDED. 

C MO  CHARGE  DENSITY  IS  COMPUTED 

C 

C XRl,  XR2,  XR3;  ARRAYS  WITH  THE  COMPONENTS  OF  THE  POSITION  VECTOR 

C ON  THE  SPHERE 

C SRF:  (SOLID  ANGLE  SUBTENDED  BY  EACH  PATCH)/(TWO  PI) 

C TF:  TIME  OF  FIRST  ARRIVAL  OF  THE  WAVE  AT  EACH  PATCH 

C CO:  COEFFICIENT  FOR  CONTRIBUTION  TO  SELF-PATCH 

C DO.  SAME  AS  CO..  FOR  THE  TIME-DERIVATIVE  TERM 

C TH2>  TH3:  COMPONENTS  OF  THE  UNIT  VECTOR  ALONG  MERIDIANS 

C PHI,  PH2:  COMPONENTS  OF  UNIT  VECTOR  ALONG  PARALLELS 

C FUl,.  FU2,  FD.1,  FD2,  FCl,  FC2:  INFORMATION  FOR  INTERPOLATION 

C THE! A:  POLAR  ANGLE  OF  PATCH 

C PHIJ;  AZIMUTHAL  ANGLE  OF  PATCH 

C DPHII;  INCREMENT  OF  AZIMUTHAL  ANGLE 

C RB.1,  RB2,  R.B3,  RB4;  ARRAYS  THAT  CONTAIN  BIT  ARRAYS 

C COJ,  COJT,  COJU,  COJTU,  COJV,  COJTV;  ARRAYS  WITH  COEFFICIENTS 

C NEEDED  TO  COMPUTE  NEIGHBORING-PATCH  CORRECTIONS 

C SINTl,  SINT2,  SI NTS:  CARTESIAN  COMPONENTS  OF  INTEGRALS 

C DRl  ..  DR2,  DR3:  COMPONENTS  OF  R = X - X ' 

C DR,  C.0RM2,  DRM3:  WORK  ARRAYS 

C CJTHD,  CJPHD:  COMPONENTS  OF  THE  RETARDED  SURFACE  CURRENT  DENSITY 

C IN  THE  THETA-  AND  PHI -DIRECT  IONS 

C DJTHD,  DJPHD:  SAME  FOR  TIME  DERIVATIVE  OF  SURFACE  CURRENT  DENSITY 

C UTH..  UPH,  AO,  Al,  A2,  A3,  A4,  YM2,  YMl,  YO,  Yl,  Y2‘  WORK  ARRAYS 

C CJTHA,  CJPHA;  COMPONENTS  OF  SURFACE  CURRENT  DENSITY  AT  TIME  T 

C CJTH,  CJPH:  TWO-DIMENSIONAL  ARRAYS  WITH  COMPONENTS  OF  ALREADY  CALCULATED 

C CURRENT  DENSITY 

C FT  INFORMATION  FOR  TIME  INTERPOLATION 

C ML,  NR,  NUI,  NU2,  NDl,  ND2;  INFORMATION  FOR  SPATIAL  INTERPOLATION 

C NET,  NPQINT;  INFORMATION  FOR  NEIGHBOR ING-PATCH  CORRECTION 

C IND:  INDEX  ARRAY  FOR  VECTOR  PROCESSING 

C IJ:  IJ(I)  I,  FOR  VECTOR  PROCESSING 

C BT.l  BIT  ARRAY  SELECTS  ALL  PATCHES  EXCEPT  THE  SELF  PATCH 

C BT2:  BIT  ARRAY  TO  INDICATE  WHETHER  THE  PULSE  HAS  REACHED  A PATCH 

C BT3,  BT4:  BIT  ARRAYS  TO  INDICATE  NECESSITY  OF  SHIFTING  INTERPOLATION 

C 

PARAMETER  (LL==LLG,  LC=LCG,  LL12=LL12G,  LL20^-LL20G. 

1  LB  ^:LL**2/64+l,  LB  1 -LL/64+1  ) 

BIT  BTl (LL), BT2(LL), BT3(LL, LL ) , BT4(LL,  LL > 

CHARACTER  TITLE^l•80,  SUBT^^SO,  XLAB*30,  YLAB*30,  QAL.-»4 
COMMON  /RCOM/  XR 1 ( LL ) , XR2 ( LL ) , XR3 ( LL ) , SRF ( LL ) , TF ( LL ) , 

1 CO(LL), DO(LL), THl (LL), TH2(LL), TH3(LL), 

2 PHI (LL), PH2(LL), FUl (LL), FU2(LL), FDl (LL), FD2(LL), FCl (2), FC2(2). 

3 THETA (LL), PHII (LL), DPHII (LL), RBI (LBl ), RB2(LB1 ), RB3(LD). RB4(LD). 

4 COJ (6, LL12) , C0JT(6, LL12), COJU (2, LL12) , C0JTU(2, LL12) - 

5 C0JV(2, LL12), COJTV (2, LL 1 2 ) , SI NTl ( LL ) , SINT2(LL), SINT3(LL). 

6 DR  1 ( LL ) , DR2  < LL ) , DR3 ( LL ) . DR ( LL ) , CDRM2 ( LL ) , DRM3 ( LL ) , 

7 CsJTHD(LL),  CJPHD  (LL),  DJTHD(LL),  DJPHD  (L.L), 

8 UTH  ( L L ) , UPH  ( LL  ) , AO  ( LL  > , A 1 < LL  ) , A2  ( LL  ) , A3  ( LL  ) , A4  ( L.L  ) , 

9 YM2 ( LL ) , YMl ( LL ) , YO ( LL ) , Yl ( LL ) , Y2 ( LL ) , 
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A CJTHA<LL),  C-JPHA(LL),  CJTH<LL,  LC ) , CJPH(LL,  LC)>  FT(LL  ), 

B NL(LL), NR(LL),  NUl (LL>,  NU2<LL)>  NDl (LL),  ND2(LL), 

C NEI  <LL,  LL20),  NPOINKLL,  LL20),  IND(LL).  IJ<LL) 

COMMON  /ICOM/  IX, IS, ISl, Nl, N2, NI, ICORR, IPLOT, INTER, IDEB, 

1 IPATCH,  N3,  NU 

COMMON  /LABEL/  TITLE, SUBT, XLAB, YLAB 

COMMON  /PARAM/  ALPHA, BETA, DT. RAD, RADMl , CDTMl , THSC, SR i , SR2, SR3 
COMMON  /CONST/  AMUO, TWOMUO,  TWOEPO,  C,  CMl,  PI,  TOPI, PITO, OS, OTW, OTF 
EQUIVALENCE  (RB1,BT1),  <RB2, BT2), (RB3, BT3),  <RB4, BT4) 

DATA  IZO/0/, BT1/LL*E ' 1 V 
WRITE (XLAB,  7) 

WRITE! YLAB,  S) 

LL2=2+^LL 
OS--=l.  /7. 

OTW--^l.  / 1 2. 

DTF:=1,  /24. 

C==3.  E8 
CM  1=1.  /C 

PIT0=2.  -itATANU.  ) 

PI-^2.  -«PITO 
TOP  1-^4,  ftp  I TO 
CDNV==PI/180. 

AMU0-'--PI-ft4.  E--7 
TW0MU0=2.  /AMUO 
TW0EP0=2.  *CM 1^*2/ AMUO 
PRINT  MINI)" 

READ  ICORR, ISKIP, IPLOT, INTER, IDEB, IPATCH. QAL 

PRINT  -ft,  '•  ICORR==^  ICORR,  S ISKIP=  S ISKIP,  ^ IPLOT^^: IPLOT, 

1 L.  INTERS  L.  INTER,  ^ IDEB=  ' , IDEB,  IPATCH=  ' , I PATCH,.  ' , QAL=  ^ QAL 

READ  IN  VALUES  OF  CONSTANTS  FOR  THE  PULSE 


READ  4,  DT, RAD, ALPHA,  BETA 
READ  5,  NI 
READ  5,  N1,N2,  N3 

IF<N1,  GT.  100)  STOP  2 
IF(W3.  EQ.  0)  N3=^5 
R ADM 1 = 1.  /RAD 
IF( ISKIP.  EQ.  1 ) 1C0RR=0 
COMPUTE  COMPONENTS  OF  VECTORS  ON 
SELF-PATCH  AND  NEIGHBORING-PATCH 
CALL  RUNT 

THEN 

I XRl ( I ) 

TF(I)  CO(I) 

TH2 ( I ) TH3 ( I ) 

PHKI)  PH2(I) 


SPHERE  AND  INFORMATION  NEEDED  FOR 
CORRECTIONS 


IF<  IDEB.  EQ.  1 ) 
WRITE (7, ft) ' 

1 'SRF(I) 

2 ^THl(I) 
WRITE <7,  ft)  ' 


XR2( I ) 
Don ) 

FUl ( I ) 


XR3( I ) 


1 ' FU2 ( I ) FD 1 ( I ) FD2 ( I ) ' 

WRITE(7, ft)'  NL(I)  NR(T)  NUl < I ) NU2(I)  NDl ( I ) ND2(I)' 
DO  90  1 = 1,  IS 

WRITE (7, 2)  I,  XRl  ( I ),  XR2( I ),  XR3( I ),  SRF < I ) , TF( I ) , 

1 C0( I >, D0( I ) , THl ( I ),  TH2< I ), 

2 T H3  ( I ),  PH  1 ( I ) , PH2  ( I ) , FU 1 < I ) , FU2  < I ) , FD  1 ( I ) .•  FD2  ( I ) 
WRITE(7,  1 ) NL(  I ),  NR(  I ),  NUl  ( I ),  NU2(  I ),  NDl  ( I ),  N.D2(  I ) 
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90 


9f-; 

C 

100 

C 

C 


CONTINUE 

IF( ICGRR.  GE,  2)  THEN 
DC)  95  IS 

WRITE  (7,  1)  I,  <NEI  < I,  J),  NEI  (I,  ,l  ) ) 

CONTINUE 
END  IF 
END  IF 

NO  CALCULATIONS  IF  NI  <'  0 
IF(NI.LT.  0)  STOP 
DO  100  I==-l,  IS 
I J( I )=I 
CONTINUE 
DTr'11--=:l.  /DT 
CDTM^=CMH«-DTM1 

COMPUTE  THE  NUMBER  OF  TIME  STEPS  TO  A DIAMETER  TO  DIMENSION  ARRAY 
NU-2,  *RAD^<-CMHtDTMl+6. 

COMPUTE  THE  NUMBER  OF  TIME  STEPS  IF  NOT  GIVEN 
IF<NI,  EQ.  0)  NI=^6.  *RAD/(C*DT) 

PRINT  'NI^MvII,  ^ N1=SN1>',  N2='.N2>'>  N3=UN3, 

1 U DT==',DT,  S RAD=-'.  RAD.  ^ ALPHA- ALPHA,  ' , DETA=',BETA 
PRINT  NU.  ' ROWS  NEEDED  OUT  OF  '.LC 
IF(NU. GT. LC)  STOP  3 
NUM2=^NU-2 

IF< ICORR.  EQ.  1 > THEN 

PRINT  Jt.  'SELF-PATCH  CORRECTIONS  ONLY' 

ELSE  IF( ICORR.  EQ.  2)  THEN 

PRINT  'NEIGHBORING-PATCH  CORRECTIONS  ONLY' 

ELSE  IF( ICORR.  EQ.  3)  THEN 

PRINT  'SELF-  AND  NEIGHBOR  ING-PATCH  CORRECTIONS' 

ELSE 

PRINT  •«■..  'MO  PATCH  CORRECTIONS' 

END  IF 

ENCODE  INFORMATION  FOR  TITLES  OF  GRAPHS 
IF( ICORR.  EQ.  1 ) THEN 
IF(  INTER.  EQ  0.^  THEN 

WRITEv'SUBT.  11  ) ALPHA.  BETA.  DT.  NI.  Nl.  N2.  N3 
ELSE 

WRITE(SUj3T,  12)  ALPHA.  BETA.  DT.  NI . Nl , N2.  N3 
END  IF 

ELSE  IF ( ICORR.  EG.  2)  THEN 
IF ( INTER.  EQ.  0)  THEN 

WRI  TEISUBT.  13)  ALPHA. BETA. DT. NI.  Nl.  N2.  N3 
ELSE 

WRITETSUBT.  14)  ALPHA.  BETA.  DT.  NI.  Nl.  N2.  N3 
END  IF 

ELSE  IFdCORR.  EQ.  3)  THEN 
IF ( INTER.  EQ.  0)  THEN 

WRrTE(SUBT.  15)  ALPHA.  BETA.  DT.  NI.  Nl.  N2.  N3 
ELSE 

WRITE(SUBT.  16)  ALPHA. BETA. DT.  NI,  Nl.  N2.  N3 
END  IF 
ELSE 

IF ( INTER.  EQ.  0)  THEN 

WRITE (SUBT. 9)  ALPHA. BETA.  DT.  NI,  Nl.  N2.  N3 
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ELSE 

WRITE(SUBT,  10)  ALPHA, BETA, DT, NI,  Nl,  M2,  N3 
END  IF 
END  IF 

CALL  Q5RQUEST(  'LFN= 'SCFL  ' //QAL,  'MXR=^  S 2400.  'LEN=^50, 

1 'RT=^  'W'  ) 

OPEN  (9,  FILE=='SCFL  V/QAL,  FORM==  ^UNFORMATTED  ' ) 

T--=--DT 

NF2^S-ttIS 

LENGTH=^N  I * I S-«-2/  512+10 
IF< ISKIP.  EG.  1 ) THEN 

C COMPUTE  FAR  FIELDS  FROM  SAVED  FILE  SCSVLQAL 

CALL  Q5ATTACH(  'LFN=^  •'SCSVV/QAD 

OPEN  (2,  FILE='SCSVV/QAL,  STATUS="OLD^  ACCESS:= 'DIRECT 
1 FORM== 'UNFORMATTED',  RECL:=NF2) 

T--NI-ttDT 
QO  TCI  400 
END  IF 

CALL  Q5RQUEST(  'LFN=--=',  'SCSV'//QAL,  ' MXR=-' ' , NF2,  'LEN= ',  LENGTH, 
1 'SFO=', 'D') 

OPEN (2,  FILE='SCSV'//QAL,  FORM-=  'UNFORMATTED', 

1 RECL=-NF2,  ACCESS=--= 'DIRECT  ' ) 

TL=SECOND( ) 

NUM=--a 
NREC^=M 
I XH----0 

/-• 

C INITIALIZE  ARRAYS 

r 

CJTHA( 1; IS)=0. 

CJPHA<  1;  IS)^^--0. 

CJTHIX  1 ; IS)-=^0, 

CJPHD(  1;  IS)=-0. 

DJTHD(  1;  IS)===0. 

DJPHDd;  IS):=0, 

UTH(1; IS)^0. 

UPH<  .1;  IS)=^0. 

DO  120  J=^1,LC 
CJTH(  1,  vJi  IS)-0. 

CJPPK  1.  Ji  IS)=0. 

120  CONTINUE 
r 

C MAIN  LOOP 

C 

DO  300  IZ--2,NI 
NUM=NUM+1 

PRINT  17,  ' OF  ',  NI 
C SHIFT  THE  NUMBERS  IN  THE  ARRAYS 

DO  170  1=2, NU 

CJTH< ] , I~li  IS)=CJTH( 1,  I;  IS) 

C JPH  < 1 , I -- 1 , IS)  =C  JPH  < 1 , I ; IS) 

170  CONTINUE 

C COMPUTE  DENSITIES  AT  EACH  PATCH  ON  THE  SPHERE 

DO  200  IX=1, IS 
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C RETRIEVE  THE  RADIUS  VECTOR  (NORMAL.) 

X=XR1 ( IX ) 

Y=-vXR2(;(X  ) 

Z~XR3( IX ) 

C FIND  THE  INCIDENT  FIELDS 

CALL  INFL(5INT10, SINT20,  SINT30, 

1 X,  Y,  Z,  T.  «-210) 

DRl  ( li  IS)r-^X-XRl  < 1;  IS) 

DR2(li  IS)=--Y-XR2(  1,  IS) 

DR3(1,  IS):--=Z-XR3(1,  IS) 

CDRM2U;  IS)'^DR1  ( .1..  IS  ) h^*2+DR2  (1 ; IS  ) **2+DR3  < 1 ; IS)**2 
DR( i;  IS)=VSQRT(CDRM2< 1,  IS); DR( li  IS) ) 

DRM3(1; IS)=CDTM1*DR< 1, IS) 

IND( 1, IS)-VINT(DRM3( 1; IS); IND( 1; IS) ) 

FT (1 ; IS) =VA I NT ( DRM3 (1 ; IS); IS ) -DRM3 ( 1 ; IS) 

C SHIFT  THE  FIVE  POINTS  FOR  INTERPOLATION  IF  NECESSARY 

WHERE<BT3( 1 , IX, IS) ) 

INDU,  IS)  = IND(  1..  IS)+1 
FT( 1;  IS)=FT( 1;  IS)  + l. 

END  WHERE 

WHERE(BT4( 1, IX; IS) ) 

IND( 1; IS)=IND< 1, IS)+2 
FT( 1 , IS)=FT< 1 ; IS) +2. 

END  WHERE 

INDd,  IS)-IJ(1;  IS)-IND(1,  IS)*LL 
,tm  ( I X ) =B  '0-' 

C DETERMINE  WHETHER  THE  PULSE  HAS  REACHED  THE  PATCH 

BT2(1;  IS)==T-CM1*DR  (1 , I S ) . LE.  TF  (1 , I S ) 

WHERE  <BT1  ( ,t  ..  IS)  ) DRM3<  1..  IS)^1.  /DR  ( 1;  IS) 

CDRM2(1;  IS)=DRM3(1;  IS)*^f2 
DRM3(1; IS)=CDRM2(1; IS ) *DRM3 ( 1 , IS) 

CDRM2(1;  IS)-CDRM2(1;  TS);«-CDTM1 
C DETERMINE  RETARDED  INTERPOLATED  FIELDS 

IF (INTER.  EQ.  O)  THEN 

CALL  INTRPOI IND, FT, CJTH, CJTHD, DJTHD, AO, Al, A2, 
t YM2, YMl, YO, Yl, Y2) 

CALL  INTRPOdND,  FT,  CJPH,  CJPHD,  DJPHD,  AO,  Al,  A2, 

1 YM2, YMl, YO, Yl, Y2) 

ELSE 

CALL  INTRPl ( IND, FT, CJTH, CJTHD, DJTHD, AO, Al, A2, A3, A4, 

1 YM2, YMl, YO, Yl. Y2) 

CALL  INTRPl  ( IND,  FT,  CJPH,  CJPHD,  DJPHD.  AO.  Al,  A2,  A3,  A4. 

.1  YM2,  YMl,  YO,  Yl , Y2) 

END  IF 

5L-T  FIELDS  TO  ZERO  IF  PULSE  HAS  NOT  REACHED  PATCH 
( INTERPOLATION  MAY  MAKE  FIELDS  NONZERO  BEFORE  THEY  ARE) 

WHERE  (BT2d,  IS)  ) 

CJTHDd;  IS)=0. 

CJPHD (1;  IS)=0. 

DJTHD d;  IS)=0. 

DJPHDd;  IS)=0. 

END  WHERE 

DETERMINE  THE  CARTESIAN  COMPONENTS  OF  THE  CONTRIBUTION  OF  THE 
PATCH  TO  THE  INTE(3RAL 
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UTH(  1,  IS)=CDRM2(  1,  I S ) «-DJTHD <1 , IS)+DRM3(1,  IS)*CJTHD(  li  IB) 
UPHd.  IS):=CDRM2(  1..  IS)#DJPHD(  li  IS)+DRM3(  li  IS)-«-CJPHD(  1;  IS) 
SINTl  ( 1;  IS)=-SRF(1;  I S ) «■  ( DR2 ( 1 , IS)4«-UTH<  1;  IS)^^TH3(  li  IS) 

1 -DR3(1; IS)*<UTH( 1 , IB)*TH2( li IS)+UPH( li IS)*PH2< 1, IS) ) ) 

SINT2(  li  IS)  = -SRF(  li  IS)-?^«:DR3(  li  IS)-«-(UTH(  li  IS)*TH1  < li  IS) 

1 +UPH(li IS)*PH1 ( li IS) )-DRl( li IS)*UTH( li IS)*TH3( 1> IS) ) 

SINT3(li  IS)=-SRF(li  IS)*(DR1(1,  IS)*(UTH(li  IS)^^TH2(li  IS) 

1 -+-UPH(li  IS)*PH2(li  IS))-DR2(li  IS)*(UTH(li  IS)-»THl(li  IS) 

2 +UPH( 1, IS)*PH1< li IS) ) ) 

SNT1=Q8SBUM(SINT1 < li IS). BTK li IS) )+SINT10 
SNT2=-Q8SSUM<BINT2<  li  IS).  BTl  < 1,  IS)  )+SINT20 
SNT3=Q3SSUM(SINT3(  li  IS).BTld.  IS))+SINT30 
BTl  (IX  )=B  d ■' 

C 

C FILL  IN  APPROXIMATE  VALUES  FOR  CURRENTS  AND  DERIVATIVES 

IF(  IX.  EQ.  1 ) THEN 
CJTHD( 1 )=SNT2 
CJPHD< 1 )=-SNTl 
CJTH( 1, NU)=SNT2 
CJPH( 1. MU)^-SNT1 
ELSE 

CJTHD(  rX)=-PHl  ( IX)«-SNT1-PH2(  IX)*SNT2 
C JPHD  a X ) =^TH1  (IX)  *SNT1  +TH2  (IX)  ^<SNT2+TH3  (IX)  *SMT3 
CJTH( IX, NU)=CJTHD( IX ) 

CJPH( IX. NU)=CJPHD( IX ) 

END  IF 

IF(  INTER.  E(3.  0)  THEN 

CALL  I NTRP2 ( 2.  . C JTH.  I X . NU-2. D JTHDX ) 

CALL  INTRP2(2. . CJPH. IX. NU-2, DJPHDX ) 

ELSE 

CALL  INTRP3(2.  , CJTH,  I X, NU-2. DJTHDX ) 

CALL  INTRP3(2.  , CJPH.  IX, NU-2. DJPHDX) 

END  IF 

DJTHD( IX JTHDX 
DJPHD( IX )=DJPHDX 

COMPUTE  CORRECTIONS  FOR  NEIGHBORING  PATCHES 

C IF  ICORR  = 0 OR  1.  NEKIX.l)  = 1 AND  LOOP  IS  SKIPPED 

DC)  .190  IXX==2.  NEI  ( IX.  1 ) 

IXP^^NEI  ( IX,  IXX  ) 

IF(T-CMt*DR( IXP).  LE.  TF( IXP)  ) GO  TO  190 
NPNTn^NPOINT(  IX.  IXX  ) 

CGJU^CQJ(  1.  NPNT) 

C0J2=C0J(2.  NPNT) 

CDJ3-“=C0J(3,  NPNT) 

CGJ4::=CGJ(4,  NPNT) 

COJ5=dlOJ(  5..  NPNT) 

C(3J6^--C0J(6,  NPNT) 

COJTi-COJTv 1,  NPNT) 

CGJT2:^CQJT(2.  NPNT) 

CO  JT3=--C0  JT  ( 3.  NPNT  ) 

CCJ  JT4  =^CO  JT  ( 4.  NPNT  ) 
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COJTS^^CQJKS.  NPNT) 

CDJT6==C0JT<6,  NPNT) 

COJUi=^COJU(  1,  NPNT) 

C0JU2=--C0JU(2.  NPNT) 

CDJTUl“-“COJTU(  1;  NPNT) 

C0JTU2-=CGJTU(2,  NPNT) 

COJVl=COJV( 1,  NPNT) 

C0JV2=CDJV(2.  NPNT) 

CDJTVl=--COJTV(  1.  NPNT) 

CGJTV2=C0JTV(2;  NPNT) 

IF(  IXP.  EQ.  .1.  OR.  IXP.  EQ.  IS)  THEN 

COMPUTE  NEIGHBORING-PATCH  CORRECTIONS  FROM  CAPS,  IF  NECESSARY 
IFdXP.EQ.  1)  THEN 
ISP  ■•=2 
ELSE 

ISPr=IS.l 
END  IF 

DETERMINE  THE  INCREMENTS  WITH  RESPECT  TO  X 

CvJlD2=^THl  ( ISP)<t-CJTHD<  ISP)+PH1  ( ISP  ) ^^C  JPHD ( I SP  ) 
CJ2D2=--TH2(  ISP  )*CJTHD(  ISP)+PH2<  ISP  ) *C  JPHD ( I SP  ) 
CJ3D2-d'H3(  ISP)-K-CJTHD(  ISP) 

N1=^--NL(  IXP  ) 

N2==NR(  IXP) 

CJlDL^dHl  (N1  )*CJTMD<N1  )+PHl  (N1  )-»-CJPHD(Nl  ) 

CJ2DL=--=TH2(N1  )*CJTHD(N1  )^-PH2(Nl  ) «-CJPHD(Nl  ) 

C-J3DL^=^TH3(N1  )-k-CJTHD(N1  ) 

CJIDR-THI (N2)«CJTHD(N2)+PH1 (N2)*CJPHD(N2) 

C J2DR==^TH2  ( N2  ) ^«-C  JTHD  ( N2  ) +PH2  ( N2  ) -«-C  JPHD  ( N2  ) 

C J3DR:--=TH3  ( N2  ) -k-C  JTHD  ( N2  ) 

IFdXP.EQ,  1)  THEN 

DXCJl=^CJlD2-(CJlDLi^FCl  ( 1 )+CJlDR*FC2(  1 ) ) 
DXCJ2=--CJ2D2-(CJ2DL*FC1  < 1 )+CJ2DR»FC2(  1 ) ) 
DXCJ3=CJ3D2-(CJ3DL^<•FC1  ( 1 ) +C J3DR*FC2  ( 1 ) ) 

ELSE 

DXCJ1=^CJ1D2-(CJ1DL*FC1  ( 2 ) +C  J1  DR^tFC2  ( 2 ) ) 
DXCv.i2=^CJ2D2~<CJ2DL-«-FCl  ( 2 ) +C  J2DR#FC2  ( 2 ) ) 
DXCJ3=CJ3D2-(C  J3DLitFCl  ( 2 ) +C  J3DR-«-FC2  ( 2 ) ) 

END  IF 

DETERMINE  THE  INCREMENTS  WITH  RESPECT  TO  X 
N1--=NU1  ( ISP) 

N2-^=NU2  d SP  ) 

CJ1PU=(TH1 (N1 )*CJTHD(N1 )+PHl ( N1 ) *C JPHD ( N1 ) )»FU1 ( ISP) 

,1.  +(TH1  (N2)*C  JTHD<N2)+PH1  (N2)#CJPHD(N2)  )<^FU2dSP  ) 

CJ2PU-^(TH2(N1  )*CJTHD(N1  )-*-PH2<Nl  )-k-CJPHD(N1  ) )*FU1  ( ISP  ) 

I + d'H2(N2)*C  JTHD  < N2 ) -<-PH2  ( N2 ) «-C  JPHD  (N2)  )*FU2(  ISP) 

CJ3PU=TH3(N1  )#CJTHD(N1  )<^FU1  ( ISP  ) 

1 +TH3<N2)^<-CJTHD(N2)*FU2dSP) 

N1=NDI ( ISP ) 

N2==ND2  (ISP) 

CJ1PD==(TH1  <N1  )^^CJTHD(N1  )+PHl  (N1  )*CJPHD(N1  ) )*FD1  ( ISP  ) 

1 +<TH1  (N2)*CJTHD(N2)+PH1  ( N2  ) *C  JPHD ( N2  ) )-«-FD2(  ISP  ) 

CJ2PD=(TH2(N1 > «C JTHD < N1 )+PH2(Nl )*C JPHD (N1 ) )#FD1 (ISP) 

1 •+■  ( TH2  < N2  ) *C  JTHD  ( N2  ) +PH2  < N2 ) *C  JPHD  ( N2  ) ) *FD2  (ISP) 

CJ3PD-TH3(N1 )#CJTHD(N1 )*FD1 ( ISP ) 
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+TH3  ( N2 ) JTHD  ( N2 ) ^^FD2  (ISP) 

DYCJl=-=CJiPU-CJlPD 
DYCJ2==CJ2PU--CJ2PD 
DYCJ3:--^CJ3PU-CJ3PD 

DETERMINE  INCREMENTS  OF  TIME  DERIVATIVES  WITH  RESPECT  TO  X 
D J 1 D2=^TH  1 ( I SP  ) *D JTHD  (ISP)  +PH 1 ( I SP  ) *D  JPHD  (ISP) 

D J2D2-^TH2  (ISP)  ^tD JTHD  (ISP)  +PH2  (ISP)  -m-D  JPHD  (ISP) 

D J3D2^^=TH3  (ISP)  -»-D  JTHD  (IBP) 

N1=--NL(  IXP) 

N2=^NR(  IXP) 

DJIDL^THI  (N1  )-J^DJTHD(Nl  )+PHl  (N1  )*DJPHD(N1  ) 
DJ2DL:--TH2(N1  ) *DJTHD  ( N1 ) +PH2  ( N1  )*DJPHD(N1  ) 

D J3DL=-^TH3  ( N 1 ) *D  JTHD  ( N 1 ) 

D J 1 DR^ TH 1 < W2 ) *D JTHD ( N2 ) +PH 1 ( N2 ) *D JPHD ( N2 ) 

D J2DR=---TH2  ( N2  ) #D JTHD  ( N2 ) +PH2  ( N2  ) *DJPHD  ( N2 ) 

D J3DR^-=TH3  ( N2  ) *D  JTHD  ( N2 ) 

IF( IXP.  EQ.  1 ) THEN 

DXDJl=--DJlD2-(DJlDL-i.^FCl  ( 1 ) +DJ1DR*FC2 ( 1 ) ) 
DXDJ2=DJ2D2-(DJ2DL#FC1 ( 1 ) +DJ2DR*FC2 ( 1 ) ) 
DXDJ3=^DJ3D2“(DJ3DL#FC1  ( 1 ) +DJ3DR*FC2 ( 1 ) ) 

ELSE 

DX D J 1 D J 1 D2  - ( D J 1 DL-t^FC  1(2)  +D J 1 DR*FC2  ( 2 ) ) 
DXDJ2-DJ2D2-\DJ2DL*FC1  ( 2 ) +DJ2DR-»^FC2  ( 2 ) ) 
DXDJ3=DJ3D2- ( DJ3DL*FC 1(2) +DJ3DR*FC2 ( 2 ) ) 

END  IF 

C DETERMINE  INCREMENTS  OF  TIME  DERIVATIVES  WITH  RESPECT  TO  Y 

N1===NU1  (ISP) 

N2=--NU2  (ISP) 

DJ1PU=(TH1  (N1  )^^DJTHD(N1  )+PHl  (N1  )*DJPHD(N1  ) )-J^FUl  ( ISP) 
1 +(TH1  (N2)^^DJTHD(N2)+PH1  (N2)*DJPHD(N2)  )-frFU2(  ISP  ) 

DJ2PU=(TH2(N1  )«DJTHD(N1  )+PH2(Nl  ) J^DJPHD(N1  ) )*FU1  (ISP) 
1 +(TH2(N2)^^DJTHD(N2)+PH2(N2)*DJPHD(N2) )*FU2( ISP ) 

DJ3PU==TH3(N1  )*DJTHD(N1  )^^FU1  ( ISP) 

1 +7  H3 ( N2 ) *D JTHD ( N2 ) *FU2 (ISP) 

N 1 ===ND  1 ( I SP  ) 

N2=^=N02  ( I SP  ) 

DJ1PD=(TH1  (N1  )-i^DJTHD(Nl  )+PHl  (N1  )*DJPHD(N1  ) )*FD1  ( ISP  ) 
;l  +(TH1  (N2)*DJTHD(N2)+PH1  ( N2 ) *DJPHD  ( N2  ) )-t^FD2(  ISP  ) 

DJ2PD=(TH2(N1  )^^DJTHD(N1  )+PH2(Nl  ) »DJPHD(N1  ) )^«FD1  (ISP) 
1 +(TH2(N2)*DJTHD(N2)+PH2(N2)*DJPHD(N2) )*FD2( ISP ) 

DJ3PD==TH3(N1  ) *DJTHD  ( N1  ) *FD1  ( ISP  ) 

1 +TH3  ( N2 ) ^^D  JTHD  ( N2 ) *FD2  (ISP) 

DYD  J 1 =--D  J 1 PU~D  J 1 PD 
D Y D J 2=^  D J2P  U - D J2P  D 

C DETERMINE  CORRECTION  TO  INTEGRAL  FROM  CAPS 

DYDJ3:-“-DJ3PU-DJ3PD 

SNT1P=^C0J3*CJPHD(  I XP  ) -C0JU2^^DXC  J3-C0JV2-K-DYC  J3 
;l  +C0JT3#DJPHD  ( I XP  ) ~CQJTU2-frDXDJ3-C0JTV2*DYD  J3 

SNT2P=-*-C0J3^fCJTHD(  I XP  ) +COJU1  *DXC  J3+C0  JV 1 *DYC  J3 
1 -C0JT3*DJTHD(  IXP)-(-CajTUl*DXDJ3  + C0JTVl*DYDJ3 

SNT3P==--C0JH^CJPHD(  IXP  )+C0J2*C  JTHD  ( I XP  ) -CO  JU  U^DXC  JPH 

1 +C0JU2#DXCJTH-C0JVl•^tDYCJPH+C0JV2^^•DYCJTH 

2 -COJTl-«-DJPHD(  IXP  )+CO  JT2»tD  JTHD  ( I XP  ) -CO  JTU 1 itDXD  JPH 

3 +C0JTU2-frDXDJTH-C0JTVm-DYDJPH+C0JTV2^H)YDJTH 
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SNT1=:SNT1-SINT1  < IXP)+SNT1P 
SNT2=SNT2-SINT2 < I XP > +SNT2P 
SNT3=SNT3“S I NT3 ( I XP ) +SNT3P 
EL.SE 

COMPUTE  NEIGHBORING-PATCH  CORRECTIONS  FROM  REGULAR  PATCHES 
DETERMINE  INCREMENT  OF  FIELDS  AND  TIME  DERIVATIVES  WITH  RESPECT  TO  THETA 
CALL  AVCUR( IXP, NUl. NU2. FUli FU2>  CJTHUO,  CJPHUO> 

1 CJTHD, CJPHD) 

CALL  AVCURdXP,  NDl,  ND2,  FDl,  FD2,  CJTHDO,  CJPHDO, 

1 CJTHDiCJPHD) 

CALL  AVCUR( IXP>  NUl,  NU2, FUl, FU2,  DJTHUO,  DJPHUO, 

;l  DJTHD.  DJPHD) 

CALL  AVCURdXP,  NDl,  ND2.  FDl,  FD2,  DJTHDO.  DJPHDO. 

.1  DJTHD,DJPHD) 

DTHCJTH=-CJTHDO-CJTHUO 
DTHCJPH=CJPHDO-CJPHUO 
DTHD  JTH==D  JTHDO-D  JTHUO 
DTHDJPH^DJPHDO-DJPHUO 
C DETERMINE  INCREMENT  OF  FIELDS  AND  TIME  DERIVATIVES  WITH  RESPECT  TO  PHI 

IXl=-NRdXP) 

IX2^WLdXP) 

CJTHXl=CJTHDdXl  ) 

CvJPHXi---CJPHDdXi  ) 

CJTHX2=CJTHDdX2) 

C JPHX2==C  JPHD  d X2 ) 

DPHC JTH=C JTHX 1 -C JTHX2 
DPHC  JPH==-C  JPHX  .1  -C  JPHX2 
DJTHXl=DJTHDdXl  ) 

DJPHXl=DJPHDdXl  ) 

DJTHX2=DJTHDdX2) 

DJPHX2=DJPHDdX2> 

DPHD JTH=D JTHX 1 -DJTHX2 
DPHD  JPH=--^DJPHX  1 -D  JPHX2 
CJTHDP=CJTHDdXP  > 

CJPHDP=-CJPHDdXP) 

DJTHDP---:D  JTHD<  IXP  ) 

DJPHDP  = DvJPHD<  IXP) 

C DETERMINE  CONTRIBUTION  TO  INTEGRAL  FROM  NEIGHBORING  PATCH 

CJITH'^C0Jl*CJTHDP+C0J4<f-CJPHDP 
DJITH:^C0JTH«-DJTHDP+C0JT4*DJPHDP 
CJIPH^=C0J2*CJTHDP+C0U5*CJPHDP 
DJIPH=COJT2^^DJTHDP+COJT5*DJPHDP 
CJIN0=C0J3*CJTHDP+C0J6-»^CJPHDP+C0JUH^DTHCJTH 
:l  +C0JU2<^DTHCJPHh-C0JV1-«-DPHCJTH<-C0JV2*DPHCJPH 

DJIN0=^^C0JT3#DJTHDP+C0JT6-t^DJPHDP+C0JTUl-«-DTHDJTH 
.1.  H-C0JTU2-S-DTHDJPH+C0JTV1-K-DPHDJTH+C0JTV2K-DPHDJPH 

CITH==-CJITH+DJITH 
CIPH^^CJIPH+DJIPH 
CINO==CJINO+DJINO 
THlP=-dHl  dXP  ) 

TH2P=  TH2 (IXP) 

TH3P:----TH3dXP  ) 

PHlP=--PH.t  < IXP  ) 

PH2P=^PH2(  IXP) 
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190 
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XP=RADMi-«XRl  < IXP) 

YP=RADM1*XR2( IXP) 

ZP=RADm-»XR3(  IXP) 

SNT 1 P=-TH  1 P-ft-C  I TH-PH 1 P #C  I PH-XP  »C  I NO 
SNT2P=-TH2P*C I TH-PH2P*C I PH-YP#C I NO 
SNT3P=-TH3P*C I TH-ZP*C I NO 
SNT1=SNT1~SINT1 < IXP)+SNT1P 
SNT2=SNT2-S I NT2 (IXP) +SNT2P 
SNT3=SNT3-S I NT3 (IXP) +SNT3P 
END  IF 
CONTINUE 

STORE  CORRECTED  VALUES  OF  INTEGRAL 
IF(  IX.  EQ.  1 ) THEN 
CJTHA( 1 )=SNT2 
CJPHA( 1 )=-SNTl 
ELSE 

C JTHA (IX) =-PH 1 ( I X > -tf-SNT 1 -PH2 (IX) »SNT2 
C JPHA  ( I X ) =TH  1 ( I X ) *SNT  1 +TH2  (IX)  -K-SNT2+TH3  (IX)  *SNT3 
END  IF 

CORRECT  FOR  SELF-PATCH  CONTRIBUTIONS  IF  REQUIRED 

IF(  ICORR.  EQ.  1.  OR.  ICORR.  EQ.  3)  THEN 

DO  BOTTOM  CAP  (THETA  = PI) 

IF(  IX.  EQ.  1 ) THEN 

C J3D2=TH3 ( 2 ) *C JTHD ( 2 ) 

N1=NL( 1 ) 

N2=NR ( 1 ) 

CJ3DL=TH3(N1 )*CJTHD(N1 ) 

C J3DR=TH3 ( N2 ) *C JTHD ( N2 ) 

DXCJ3=CJ3D2-(CJ3DL*FC1 ( 1 )+CJ3DR*FC2( 1 ) ) 

N1=NU1 ( 1 ) 

N2==NU2  ( 1 ) 

CJ3PU=TH3(N1 )*CJTHD(N1 )*FU1 ( 1 )+TH3(N2)*CJTHD(N2)*FU2( 1 ) 
N1=ND1 ( 1 ) 

N2=ND2 ( 1 ) 

CJ3PD=TH3(N1 )*CJTHD(N1 )*FD1 ( 1 ) +TH3 ( N2 ) *C JTHD ( N2 ) »FD2 ( 1 ) 
DYCJ3=CJ3PU-CJ3PD 
D J3D2=TH3 ( 2 ) *D JTHD ( 2 ) 

N1=NL( 1 ) 

N2=NR ( 1 ) 

D J3DL=TH3 ( N 1 ) *D JTHD ( N 1 ) 

DJ3DR=TH3 ( N2 ) *D JTHD ( N2 ) 

DXDJ3=DJ3D2-(DJ3DL*FC1 ( 1 ) +DJ3DR*FC2 ( 1 ) ) 

N1=NU1 ( 1 ) 

N2=NU2( 1 ) 

D J3PU=TH3 ( N 1 ) *D JTHD ( N 1 ) «FU 1(1) +TH3 ( N2 ) *D JTHD ( N2 ) »FU2 ( 1 ) 
N1=ND1 ( 1 ) 

N2=ND2( 1 ) 

D J3PD=TH3 ( N 1 ) *D JTHD ( N 1 ) *FD 1(1) +TH3 ( N2 ) *D JTHD ( N2 ) *FD2 ( 1 ) 

DYDJ3=DJ3PU-DJ3PD 

PC=CO( 1 ) 
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PD^-=DO(  1 ) 

PDD=DO(IS) 

CJTHA(  1 )=CJTHAa  > ^«-PC+DJTHD  ( 1 ) 1 25*DXC J3-PDD^^DXDJ3 

CJPHA(  1 )=CJPHA<  1 )<^PC+DJPHD(  1 )^tPD-.  1 25^^DYC  J3-PDD^f•DYDJ3 
ELSE  IF (IX.  NE,  IS)  THEN 
C 

C DO  F^EGULAR  PATCF^ES 

C 

C JTHA  < I X > =C JTHA ( I X ) * ( 1 . -CO (IX)) -DJTHD (IX) «DO (IX) 

CJPHA(  IX  )=CJPHA(  IX)  «•(  i.  +CO(IX)  ) +DJPHD ( I X ) *DO ( I X ) 

EI.SE 

P 

C DO  TOP  CAP  IF  RE(3UIRED 

C 

CJ3D2=TH3( ISl )*CJTHD( ISl ) 

N1=NL( IS) 

N2=^NR  (IS) 

C J3DL=^TH3  ( N1  ) »C  JTHD  ( N 1 ) 

C J3DR=TH3 ( M2 ) *C JTHD ( N2 ) 

DXCJ3=^-CJ3D2-(CJ3DL»-FC1  (2)+CJ3DR*FC2(2)  ) 

N1==NU1  ( IS) 

N2=--NU2  (IS) 

CJ3PU=^TH3(N1  )-«-CJTHD(Nl  ) »-FUl  ( IS  ) +TH3  ( N2  ) *C  JTHD  ( N2  ) *FU2  ( IS) 
N1=^ND1  ( IS) 

N2=ND2( IS) 

CJ3PD:“-=TH3(N1  )^^CJTHD(N1  )*FD1  ( I S ) +TH3  ( N2  ) *C  JTHD  ( N2  ) *FD2  ( IS) 

DYCJ3=^CvF3PU-CJ3PD 

DJ3D2'--^TH3(  ISl  )*DJTFHD(  ISl  ) 

N1=NL( IS) 

N2==NR  (IS) 

DJ3DL:^TH3(N1  )-«-DJTHD(Nl  ) 

D J3DR---TH3  ( N2 ) JTHD  ( N2 ) 

DXD J3-- DJ3D2-  ( DJ3DL*FC  1(2)  +DJ3DR-»^FC2  ( 2 ) ) 

N1==NU1  ( IS) 

N2^==NU2  (IS) 

DJ3PU^TH3(N1  )^^DJTHD(N1  )-»-FU,t  ( IS  ) + TH3  ( N2  ) *DJTHD  ( N2  ) -B-FU2  ( IS) 
Nl^-NDl  (IS) 

N2^=ND2(  IS) 

DJ3PD:--=TH3(N1  )^^-DJTHD(Nl  )-«-FDl  ( IS  ) +TH3  ( N2  ) *DJTHD  ( N2  ) *FD2  ( IS) 
DYDJ3==DJ3PU-DJ3PD 

CJTHA(  IS)=CJTHA(  I S ) -tt-PC+D  JTHD  ( IS)^«-PD+.  1 25*DXC  J3+PDDi^DXD  J3 
CJPHA(  IS)=CJPHA(  IS)*PC4-DJPHD(  IS)^^PD+.  1 25*DYC J3+PDD^^DYD J3 
END  IF 
END  IF 
200  CONTINUE 

210  IXH=IX-1 

C STORE  NEW  VALUES  OF  FIELDS 

CJTH( 1, NUi  IS)=CJTHA( 1;  IS) 

CJPHC 1, NU; IS)=CJPHA( 1; IS) 

IF(NUM.  EQ.  NU)  THEN 

C WRITE  RESULTS  OUT  TO  DISK  IF  ARRAY  IS  ALL  NEW 

DO  2S0  1=1, NU 

WRIT£(2, REC=NREC ) CJTH(l,IiIS) 

WRITE(2,  REC=NREC  + 1 ) CJPHd,  li  IS) 
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NREC^NREC+2 
280  CONTINUE 

NUM=^0 
END  IF 

C INCREMENT  TIME 

T=T+DT 
300  CONTINUE 

TL=SECOND( )-TL  MAIN 

PRINT  'EXECUTION  TIME  FOR  THE  LOOP  IS  ',TL,  ' SECONDS' 
C WRITE  REMAINING  RESULTS  OUT  TO  DISK 

DO  310  I=NU-NUM+1, NU 

WR  J.TE<2,  REC-=NREC)  CJTH<1,I;IS) 

WRITE (2, REC=NREC+1 ) CJPH( 1, li IS) 

NREC=NREC+2 
310  CONTINUE 

C MAKE  DISK  FILE  PERMANENT 

CLOSE < 2) 

CALL  Q5PURGE(  'LFN=- 'SC3V'//QAL,  'STATUS=  ' , ISTAT ) 
IFdSTAT.  GT.  0.  AND.  ISTAT.  NE.  1402)  STOP  4 
CALL  Q5DEFINE(  'LFN=',  'SCSV'//QAD 

OPEN  (2,  FILE='SC3V'//QAL,  STATUS=- 'OLD  ' , ACCESS^ 'DIRECT  ' , 

1  FORM^ ' UNFORMATTED ' , RECL=NF2 ) 

400  READ(*,  4,  END^-500)  THSC>  PHSC 
ITH=THSC 
lPHI:--:-PHSC 

WRITE(TITLE, 6)  ITH,  IPHI 

THSC=THSC*CONV 

PHSC=PHSC^«-CONV 

C COMPUTE  COMPONENTS  OF  SCALED  VECTOR  TO  FIELD  POINT 

SR  1 :=^S  I N ( THSC  ) -J^COS  ( PHSC  ) *CDTM  1 
SR2=^S  I N ( THSC  ) -J^S  I N ( PHSC  ) *CDTM  1 
SR3==C0S(THSC  )*CDTM1 
C COMPUTE  FAR  FIELDS 

CALL  FARFL 
GO  TO  400 
500  CL0SE(9) 

C SAVE  FILE  WITH  COMPUTED  FIELD  INTENSITIES 

CALL  G5PURGE<  'LFN=',  'SCFL'//QAL,  'STATUS=  ISTAT  ) 
IFdSTAT.  GT.  0.  AND.  ISTAT.  NE.  1402)  STOP  5 
CALL  Q5DEFINE<  'LFN=- ',  'SCFL'//QAD 
STOP 

1 FORMAT < 20 ( IX,  15)  ) 

2 FORMATS  IX.  14, 3X,  10(2X,  IPEIO.  3) ) 

4 FORMAT (El  2.  5) 

5 FORMAT (314) 

6 FORMAT ( 'PULSE  SCATTERED  IN  THE  DIRECTION  THETA= 

1 ',  PHI=  ',  15) 

7 FORMAT  ( 'TIME  i«-C  IN  METERS') 

3 FORMAT (' INTENSITY  ') 

9 r np|viAT(  'LST.  SQ,  ALPHA=^  ' . F4.  1 , ' , BETA=  ' , F4.  1 , ' , 

i WT--=',I3,  '.  N1=',I2,  ',  N2=^',I2,  N3=',I2,  ', 

10  EORMAT(  'INTERP.  , ALPHA= ' , F4.  1 , ' , PETA~  ' , F4.  1 , ' , 

1 NT=^',  13,  ',  Nl=-  ',  12,  ',  N2-^-'  ..  12,  '.  N3=  ' . 12, 

11  FGRMAK  LST.  SQ,.  ALPHA^^  ' , F4.  1 , ' , BETA^-- ' , F4.  1 , ' , 


',  14, 


DT=  ',  1PE7  1, 
NO  CORR ' ) 
DT=  ',  1PE7.  1, 
NO  CORR ' ) 
DT-',  1PE7.  1. 


12 


13 


14 

15 

16 


C 


C 

c 


c 

c 


c 


1 13,  ^ Nl-»',  12,  S N2==S  12.  ^ N3=  ^ 12,  ' 

FORMAK  'INTERP.  , AL.PHA=  F4.  1,  ^ BETA=^',F4.  1,  ' 
1 NT==',I3,  Nl^',  12,  N2-^',I2,  N3=',I2,  ' 

FORMAK'LST.  SQ,  ALPHA™ F4.  1 , BETA=^F4.  1,' 
1 NT==',.T3,  Nl==^12,  N2™',12,  N3=:-SI2,  '' 

FORMAT  ( 'INTERP.  , ALPHA™',  F4.  1,  ',  BETA=^',F4.  1,  ' 
1 NT=-^',  13,  ',  Ni=',I2,  ',  N2=-',I2,  ',  N3=--=',I2,  ' 

rORMAT<  'LST.  SQ,  ALPHA=  ' , F4.  1 , ' , BETA'^ ' , F4.  1 , ' 
1 ',  N1=',I2,  ',  N2=^',I2,  ',  N3=',.T2,  ' 

FORMAT ( 'INTERP.  , ALPHA===SF4.  1,  ',  DETA=  ',F4.  1,  ' 
i NT=--',  13,  ',  Nl=',  12,  ',  N2=^',  12,  ',  N3=',  .12,  ' 

f:md 


SLF  PCH'; 
DT*',  1PE7.  1, 
SLF  PCH') 
DT-',  1PE7.  1, 
NEI  PCH') 
DT«',  1PE7.  1, 
ME I PCH') 
DT*',  1PE7.  1, 
SLF  S<  N') 
DT=',  1PE7.  1, 


SLF  ?/  N') 


SUBROUTINE  AVCUR  ( IN,  NXl , NX2,  FXl , FX2,  CJTHfU  C-..fFHP, 
1 CJTHX.CJPHX) 

TO  COMPUTE  CURRENTS  AT  DESIRED  PHI  BY  AVERAGING 
PARAMETER  ( LL=LLG,  LC^^LCG  ) 

DIMENSION  NXl (LL), NX2(LL), FXl (LL), FX2(LL), 

1 CJTHX(LL), CJPHX(LL) 

M1==NX1  ( IN) 

N2:=NX2(  IN) 

F1=FX1 (IN) 

F2-=^F  X2  (IN) 

IF(N1.EQ.  N2)  THEN 
IF(N1.EQ.  1)  THEN 
FOR  BOTTOM  CAP 

CJTHP  =-F1  «-CJTHX  (N1  > ~F2^5-C  JPHX  ( N 1 ) 

ELSE 

FOR  TOP  CAP 

CJTHP=--F1*CJTHX  (N1  )+F2-K-CJPHX(Nl  ) 

END  IF 

CJPHP==-F2#CJTHX  <N1  ) +F 1 *C  JPHX  ( N 1 ) 

ELSE 

FDR  OTHER  PATCHES 

CJTHP  =^-F1  »CJTHX  (N1  ) +F2*C  JTHX  ( N2  ) 

CUPHP=Fl-«-C  JPHX  (N1  )+F2-K-CJPHX  (N2> 

END  IF 
RETURN 
END 

SUBROUIINE  FARFL 


C THIS  SUBROUTINE  COMPUTES  THE  RADIATION  FIELDS  AT  TIMES  IN  T1 

C IN  THE  DIRECTION  GIVEN  BY  SRI,  SR2,  SR3. 

C TMP:  TIME  ARRAY 

C AMP:  AMPLITUDE  ARRAY 

r 

PARAMETER  (LL=LLO,  LC=LCG,  LL12=^LL12G,  LL20=LL20G, 

1 LB==LLit*2/64+l,  LB  l=LL/64-«-l  ) 

D I MENS I ON  TMP ( 300 ) , AMP ( 300 ) 

COMMON  /RCOM/  XR 1 < LL  ) , XR2 ( LL > , XR3 < LL ) , SRF ( LL ) , TF < LL ) . 

1 CO(LL),  DO(LL),  THl (LL), TH2(LL), TH3(LL), 

2 PHI (LL), PH2(LL), FUl (LL), FU2(LL), FDl (LL), FD2(LL). FCl (2), FC2(2), 
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3 THETA (LL),  PHI  I (LL) , DPHI I< LL ) , RB 1 < LB  1 ) , RB2<LB1 ) > RB3(LB) , RB4(LB ), 

4 C0J<6,  LL12),  C0JT(6,  LL12),  C0vJU(2,  LL12),  C0JTU(2,  LL12), 

5 C0JV(2,  LL12),  C0JTV(2.  LL12),  SINTKLL),  SINT2(LL),  SINT3(LL)> 

6 DR  1 ( LL) . DR2  < LL) , DR3 ( LL) , DR (LL) . CDRM2  < LL) , DRM3 ( LL)  , 

7 C JTHD ( LL) , C JPHD ( LL) > DJTHD ( LL) , DJPHD ( LL) > 

8 UTH(LL),  UPHELD,  AO  < LL)  > A1  ( LL) , A2  < LL) , A3  ( LL) , A4  ( LL) , 

9 YM2<LL), YMl (LL),  YO<LL),  Y1 (LL ),  Y2<LL), 

A C JTHA ( LL ) , C JPHA ( LL ) , C JTH ( LL,  LC ) , C JPH ( LL,  LC ) , FT ( LL ) , 

B NL(LL), NR(LL), NUl (LL), NU2(LL),  NDl (LL),  ND2(LL), 

C NEI (LL, LL20), NPGINT(LL,  LL20) , IND(LL),  IJ(LL) 

COMMON  /ICOM/  IX, IS, ISl, Nl, N2, NI, ICORR, IPLOT, INTER, IDEB, 

1 IPATCH, N3, NU 

CHARACTER  TITLE*80, SUBT*80, XLAB*30, YLAB*30 
COMMON  /LABEL/  TITLE, SUBT, XLAB, YLAB 

COMMON  /PARAM/  ALPHA, BETA, DT, RAD, RADMl,  CDTMl,  THSC,  SRI,  SR2,  SR3 
COMMON  /CONST/  AMUO, TWOMUO, TWOEPO, C,  CMl,  PI,  TOPI,  PITO, OS,  OTW,  OTF 
PRINT  *,  'START  FAR-FIELD  COMPUTATION' 

CDTM2=CDTMl*i^-2 
RCDTMi=--RAD^J-CDTMl 
FAC=^(1.  E-7*T0PI*RADM1  )**2 
NU=1 

NI2=NI*2 
READ(-»,  1 ) Tl,  T2 

COMPUTE  FIRST  TIME  OF  ARRIVAL  IF  Tl  = 0 

IF(T1.  EQ.  0.  ) T1=RAD*CM1*( 1.  -2.  *SIN(.  5*THSC ) ) -DT 

COMPUTE  LAST  TIME  IF  T2  = O 

IF(T2.  EQ.  0.  ) T2=-(NI -5)*DT-RAD#CM1 

T1DT--^T1/DT 

READ(*,2)  N 

IF(N.  GT.  300)  STOP  6 

AMP(  i;  300)=— 1. 

DT1  = (T2-T1  ) / ( (M-1  )-«-DT) 

NMAX---T 1 DT+RCDTM 1 +4. 

NMIN=T lDT-RCDTMi-3. 

NTOT=^NMAX-NMIN+i 
IF(NTOT,  GT.  LC ) THEN 

PRINT  'NTOT=  ',  NTOT,  ' GREATER  THAN  LC~  ' . LC 
STOP  7 
END  IF 

IF(NMIN.  LE.  0)  THEN 
DO  130  Jl=^l,  -NMIN+1 
CJTH(  1,  Jli  IS)=0. 

CJPH( 1, Jl;  IS)^0. 

.130  CONTINUE 
WREC-0 
ELSE 

NREC=(NMIN-1 )*2 
Jl=^l 
END  IF 

DO  150  J=J1,NT0T 
NREC=NREC+2 
IF(NREC.  GT.  NI2)  THEN 

PRINT  «•,  'TIME  LIMIT  ON  SPHERE  EXCEEDED  DURING  INTEGRATION' 
PRINT  'T2=  ',T2,  ',  T==  ',T1DT*DT,  '.  K=-  ',K,  ' OUT  OF  ',N 
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STOP  S 
END  IF 

READ(2,  REC=NREC-1  ) CJTHd,  J;  IS) 

READ(2,  REC-NREC)  CJPHd,  J;  IS) 

150  CONTINUE 

DO  200  K=1,N 

NX==TlDT+RCDTm4-4. 

NN~TlDT-RCDTMl-3. 

IF(NX~NN+1.  GT.  LC)  THEN 

PRINT  'NEW  RANGE  SNX-NN+l,  ' GREATER  THAN  LC=',LC 
PRINT  *,  'NN=',NN,  S NX^=^SNX,  K-',K 

STOP  9 
END  IF 

IF<NX-NMIN4-1.  GT.  LX  ) THEN 
NA^NN-NMIN 
DO  160  L-1,  NMAX-NN+1 
MA=NA+1 

CJTH( 1, L; IS)=CJTH< 1, NA; IS) 

CJPHd.  L>  IS)=CJPHd,  NA;  IS) 

160  CONTINUE 

NNIN==NN 
ELSE 

L^NMAX-NNINf-2 
END  IF 

IF(NX.  GT.  NMAX  ) THEN 
DO  170  J==^L,  LC 
NREC=NREC+2 
IF(NREC.  GT.  NI2)  THEN 

PRINT  «->  'TIME  LIMIT  ON  SPHERE  EXCEEDED  DURING  INTEGRATION' 
PRINT  'T2=  ',T2,  T=  ',T1DT»-DT,  K=  '.K.  ' OUT  OF  ',N 
STOP  10 
END  IF 

READ (2, REC=NREC-1 ) CJTH< 1, J; IS) 

READ(2,  REC=NREC)  CJPHd,J;IS) 

170  CONTINUE 

MMAX=^-TviMIN4-L.C~l 
END  IF 

NMINl^NMIN-1 
TMP  (K)--T1DT^^DT»^C 

DR  d;  IS)=-T1DT  + XR1  ( 1;  I S ) »^SR  1 -*-XR2  < li  IS  ) -»-SR2+XR3  ( 1 ; IS)»SR3 
I ND  ( 1 I S ) I NT  ( DR  d i IS);  I ND  d ; IS)) 

FTd;  IS)=DR  d>  IS)  -VAINT(DR  < 1;  IS);  IS) 

INDd;  IS)  = I Jd,  IS)  + dNDd;  IS)  -NMINl  )*LL 
IF (INTER.  EQ.  0)  THEN 

CALL.  INTRPO(  IND,  FT,  C JTH,  CJTHD,  DJTHD,  AO,  Al,  A2. 

1 YM2, YMl, YO, Yl, Y2) 

CALL  INTRPO( IND,  FT, CJPH, CJPHD, DJPHD. AO, Al, A2, 

1 YM2,  YMl,  YO,  Yl,  Y2) 

ELSE 

CALL  INTRPl ( IND,  FT,  CJTH, CJTHD,  DJTHD,  AO, Al, A2, A3, A4, 

1 YM2, YMl, YO, Yl, Y2) 

CALL  INTRPl ( IND,  FT, CJPH, CJPHD, DJPHD, AO, Al, A2, A3, A4, 

YM2,  YMl,  YO,  Yl, Y2) 

END  IF 
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C1=Q8SSUM(  <DJTHDU.;  IS)*TH1  (1;  IS)+DJPHD(li  IS)^^PH1  (1;  IS)  ) 

1 *SRF< li IB) ) 

C2=^QSSSUM(  <DJTHD(  1;  IS)^^TH2(1;  IS)+DJPHD(1;  IS)»-PH2(li  IS)  ) 

1 ^tSRF<  1;  IS)  ) 

C3-^QSSSUM(DJTHD(  1;  IS)*TH3U;  IS)*SRF<  1;  IS)  ) 

AMP  ( K ) = ( ( C 1 *^»•2+C2■fr*2+C3^^■*2 ) i^CDTM2-  ( C 1 *SR  1 +C2-JtSR2 
1 +C3^^SR3)**2)*FAC 

T1DT=T1DT+DT1 
200  CONTINUE 

WRITE (9)  IMP 
WRITE < 9)  AMP 

IFdPLOT.  EQ.  1.  OR.  IPLOT.  EQ  2) 

1 CALL  DRAWl (1,  16,  12,  56, 00, N, 0.  , 0.  , XLAB,  YLAB. TITLE, SUBT, TMP, AMP ) 
IF(  IPLOT.  EQ.  0.  OR.  IPLOT.  EQ.  2) 

1 CALL  DRAW2<  1,  16,  12,  56,  SO,  N,  0.  , 0.  , XLAB,  YLAB,  TITLE,  SUBT,  TMP,  AMP) 
RETURN 

1 FORMAT < El 2.  5) 

2 FORMAT (21 4) 

END 

C 

SUBROUTINE  INFL(SINT1, SINT2, SINT3,  X,  Y,  2, T, *) 

C 

C THIS  SUBROUTINE  PROVIDES  THE  INCIDENT  FIELDS 

r 

COMMON  /PARAM/  ALPHA,  BETA,  DT,  RAD,  RADMl,  CDTMl,  THSC,  SRI, SR2, SR3 

COMMON  /CONST./  AMUO,  TWOMUO,  TWOEPO,  C,  CMl,  PI,  TOPI,  PI  TO,  OS,  OTW,  OTF 

2CT=^2-C^fT+RAD 

IF(  2CT.  GE.  0.  ) RETURN  1 

AZ:=ALPHA^<-2CT 

BZ=BETA*ZCT 

IF(BZ.  GE,  -25.  ) THEN 

Bl-^-BZ*(EXP(AZ)-EXP(BZ)  ) 

E2=^-B1#C 

ELSE 

BL“^0. 

£2^0. 

END  IF 

SINT1=TW0MU0KD1 

SINT2==0. 

SI NT 3=0. 

RETURN 

END 

C 

SUBROUTINE  INTRP0< IND, FT, ARR, CA, DA, AO, Al, A2, 
i YM2,  YMl,  YO,  Yl,  Y2) 


LEAST  SQUARES  INTERPOLATION 
PARAMETER ( LL=LLG, LC=LCG ) 

COMMON  /ICOM/  IX, IS, ISl, Nl, N2, NI, ICORR, IPLOT, INTER, IDEB, 

1 IPATCH, N3, NU 

COMMON  /CONST/  AMUO, TWOMUO, TWOEPO, C, CMl, PI, TOPI, PITO, OS, OTW. OTF 
DIMENSION  ARR(LL,  LC ) , CA(LL), DA(LL), YM2(LL), YMl <LL) , YO(LL), 

1 Y 1 ( LL ) , Y2 ( LL ) , AO ( LL ) , A 1 ( LL ) , A2 ( LL ) , I ND ( LL ) , FT ( LL ) 
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Y2( 1;  IS)=Q8VGATHR(ARR(1,  NU+2i  1 ),  IND< 1;  IS); Y2< 1;  IS) ) 

Y1  ( 1;  IS)==QSVGATHR(ARR(  1,  NU+l;  1 ),  IND<  1;  IS);  Y1  (1;  IS)  ) 

Y0(  1;  IS)=QSVGATHR(ARR(  1,  NU;  1 ),  INDd;  IS);  Y0<  1;  IS)  ) 

YMl  ( 1;  IS)=-~Q8VGATHR(ARR<  1,  NU-1;  1 ),  INDd;  IS);  YMl  d;  IS)  ) 

YM2d;  IS)-Q8VGATHR<ARRd,  NU-2;  1 )/  INDd;  IS);  YM2d;  IS)  ) 

AOd;  IS)  = (YM2d;  IS)+YM1  d;  IS)+YOd;  IS)+Y1  d;  IS)+Y2d;  IS)  )*.  2 
Aid;  IS)-(Y2d;  IS)+.  5^^(Y1  d;  IS) -YMl  d;  IS)  )-YM2d;  IS)  )*.  2 
A2d;  IS)  = (Y2d;  IS)-.  5*(Y1  d;  IS)+YM1  d;  IS)  )-YOd;  IS)+YM2d;  IS)  )*OS 
CAd;  IS)=AOd;  IS)+FTd;  IS)^^A1<1;  IS) 

1 +<FTd;  IS)*-K^2-2.  )*A2d;  IS) 

DAd;  IS)*A1  <1;  IS)+FTd;  IS) *2.  *A2d;  IS) 

RETURN 

END 

C 

SUBROUTINE  INTRPl ( IND, FT, ARR, CA,  DA,  AO,  Al,  A2,  A3,  A4, 

1 YM2, YMl. YO,  Yl,  Y2) 

C 

C POLYNOMIAL  INTERPOLATION 

C 

PARAMETER  < LL’^LLG,  LC=LCG  > 

COMMON  /ICOM/  IX,  IS,  ISl,  Nl,  N2,  NI..  ICORR,  IPLOT,  INTER,  IDEB, 

1 IPATCH, N3,  NU 

COMMON  /CONST/  AMUO,  TWOMUO,  TWOEPO,  C,  CMl,  PI,  TOPI,  PI TO,  OS, OTW,  OTF 
DIMENSION  ARR(LL, LC), CA<LL), DA(LL>, YM2(LL),  YMl (LL),  YO(LL), 

1 Yl  (LL),  Y2<LL),  AO<LL),  AKLL),  A2(LL),  A3(LL).  A4(LL), 

2 IND(LL),  FT(L.L) 

Y2d;  IS)=Q8VGATHR(ARRd,  NU+2;  1 ),  INDd;  IS);  Y2d;  IS)  ) 

Yl  d;  IS)==Q8VGATHR(ARRd,  NU+l;  1 ),  INDd;  IS);  Yl  d;  IS)  ) 

YO(  1;  IS)=QSVGATHR(ARRd,  NU;  1 ).  INDd;  IS);  YOd;  IS)  ) 

YMl  d;  IS)==Q8VGATHR(ARR(1,  NU-1;  1 ),  INDd;  IS);  YMl  d;  IS)  ) 

YM2d;  IS)==Q8VGATHR(ARRd,  NU-2;  1 ),  INDd;  IS);  YM2d;  IS)  ) 

C AOd;  IS)=YO(  1;  IS) 

Aid;  IS)=OTWK-(  YM2d;  IS)-Y2d;  IS) -8  * ( YMl  d ; IS ) -Yl  (1 ; I S ) ) ) 

A2d;  IS)=0TF*(-YM2d;  IS)~Y2d;  IS) + 16.  ■«■  ( YMl  d ; I S ) +Y 1 (1 ; IS)  ) 

1 -30.  -^fYOdi  IS)  ) 

A3(  1,  IS)=0TW-K-(-YM2d;  IS)+Y2(  1;  IS) +2.  *(  YMl  ( 1;  IS)-Y1  d,  IS)  ) ) 

A4(  1;  IS)=OTF*(  YM2(  1;  IS)+Y2d;  IS) -4.  *(YM1  d;  IS)+Y1  ( 1;  IS)  ) 
i +6.  *Y0( 1; IS) ) 

CAd,  IS)=YOd;  IS)+FTd;  IS)*(Ald;  IS)+FTd;  IS) 

1 ■»■  ( A2  d ; I S ) +FT  d ; I S ) «■  ( A3  d ; IS)  +FT  (1 ; IS)  -«-A4  d ; IS)))) 

DAd;  IS)==A1  d;  IS)+FTd;  IS)*(2.  *A2d;  IS)+FT(1;  IS) 

1 ^(3.  *A3d;  IS)+FTd;  IS) *4.  #A4d;  IS)  ) ) 

RETURN 

END 

SUBROUTINE  INTRP2(TT, ARR, I, NTP, DA) 

C 

C LEAST-SQUARES  INTERPOLATION,  SCALAR  VERSION,  DERIVATIVE  ONLY 

r 

PARAMETER  (LL  ^^LLG,  LC=LCG) 

DIMENSION  ARR(LL, LC) 

COMMON  /CONST/  AMUO, TWOMUO, TWOEPO,  C,  CMl,  PI,  TOPI,  PI  TO,  OS.  OTW.  OTF 
YM2=ARRd,  NTP-2) 

YM1=--^ARR  (I,  NTP-1  ) 

YO^ARRd,  NTP) 
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Y1=ARR( I>  NTP+i ) 

Y2=ARR( I. NTP+2) 

C A0==  ( YM2+YM 1 +YO+Y 1 +Y2 ) *.  2 

A1  = (Y2+.  5->^(Yl-YMl  )-YM2)*.  2 
A2=<Y2-.  5*<Y14-YM1  )-Y0+YM2)*0S 
C CA=A0+TT*Al  + (TT*s^2-2.  )*A2 

DA=AH-TH^2.  *A2 
RETURN 
END 
C 

SUBROUTINE  INTRP3(TT, ARR, I, NTP, DA) 

POLYNOMIAL  INTERPOLATION,  SCALAR  VERSION,  DERIVATIVE  ONLY 

PARAMETER  ( LL==LLG,  LC=LCG ) 

DIMENSION  ARR(LL,  LC) 

COMMON  /CONST/  AMUO. TWOMUO,  TWOEPO, C, CMl , P I , TOP  I , P I TO, OS,  OTW, DTP 
YM2=ARR( I,  NTP-2) 

YM1=ARR( I, NTP-1 ) 

YO=ARR< I,  NTP) 

Y.l=ARR(  I,  NTP+1  ) 

Y2=ARR( I,  NTP+2) 

Al=OTW*( YM2-Y2-8.  -fr<YMl-Yl ) ) 

A2=0TF*(-YM2-Y2+16.  YMl+Yl ) -30.  *Y0) 

A3=0TW^^(-YM2+Y2+2.  *(YM1-Y1  ) ) 

A4=0TF-k-(  YM2+Y2-4.  *(YMl+Yl)+6.  -«-Y0) 

C CA=YO+FT*< A1+TT*(A2+TT*(A3+TT#A4) ) ) 

DA=Al-i-TT<^(2.  ->tA2-»-TT^t<3.  *A3+TT^^4.  *A4 ) ) 

RETURN 

END 

SUBROUTINE  RUNT 
C 

C THIS  SUBROUTINE  COMPUTES  COMPONENTS  OF  THE  RADIUS  VECTOR 

C AND  FACTORS  NEEDED  FOR  CORRECTIONS 

C 

C N1 -NUMBER  OF  STRIPS  ON  THE  SPHERE  (VALUES  OF  THETA) 

C N2:=MAXIMUM  NUMBER  OF  PATCHES  IN  A STRIP  (VALUES  OF  PHI  ON  EQUATOR) 

C N3=MINIMUM  NUMBER  OF  PATCHES  IN  A STRIP  (NEXT  TO  A POLAR  CAP) 

C NPHI  IS  THE  ARRAY  WHERE  THE  NUMBER  OF  PATCHES  IN  A STRIP  IS  STORED 

C IS=TOTAL  NUMBER  OF  ANGLE  PAIRS 

C 

PARAMETER (LL=LLG, LC=LCG, LL12=LL120,  LL20=LL20G, 

1  LB=^LL**2/64+l,  LBl=LL/64+l  ) 

DIMENSION  NPHI ( 100) 

BIT  BT3(LL, LL), BT4(LL,  LL) 

EXTERNAL  F003,  F103.  F013.  F203,  FU3,  F023 
EXTERNAL  F002, F102, F012,  F202,  FI  12,  F022 
COMMON  /RCOM/  XR 1 ( LL ) , XR2 ( LL ) , XR3 ( LL ) , SRF ( LL ) , TF ( LL ) , 

1 CO(LL), DO(LL). THl (LL), 7H2(LL),  TH3(LL), 

2 PHI (LL), PH2(LL),  FUl (LL),  FU2(LL),  FDl (LL),  FD2(LL),  FCl (2),  FC2(2), 

3 THETA (LL), PHII (LL),  DPHII (LL),  RBI (LBl ),  RB2(LB1 ),  RB3(LB),  RB4(LB), 

4 C0J(6, LL12), C0JT(6, LL12),  C0JU(2,  LL12),  C0JTU(2,  LL12), 

5 C0JV(2, LL12), C0JTV(2,  LL12 ) , SINTl ( LL ) , SINT2(LL), SINT3(LL), 

6 DR  1 ( LL ) , DR2 ( LL ) , DR3 ( LL ) , DR ( LL ) , CDRM2 ( LL ) , DRM3 ( LL ) , 
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7 C JTHD ( LL) , C JPHD < LL) , D JTHD < LL) , D JPHD < LL) . 

3 UTH<LL), UPH(LL), AO ( LL ) , A1 < LL) , A2 ( LL) , A3 ( LL) , A4 < LL ) . 

9 YM2<LL), YMl <LL)/ YO ( LL ) , Yl( LL ) > Y2 ( LL ) , 

A C JTHA  < LL ) , C JPHA ( LL ) > C JTH  < LL,  LC ) , C JPH  < LL, LC ) , FT ( LL ) , 

B NL(LL), NR(LL), NUl (LL), NU2(LL), NDl (LL), ND2(LL), 

C NEI  (LL,  LL20),  NPOINKLL,  LL20),  IND(LL),  IJ(LL) 

COMMON  /ICOM/  IX, IS, ISl, Nl, N2, NI, ICORR, IPLOT, INTER, IDEB, 

1 IPATCH, N3, NU 

COMMON  /PARAM/  ALPHA, BETA, DT, RAD, RADMl , CDTMl , THSC,  SR  1 , SR2,  SR3 
COMMON  /CONST/  AMUO, TWOMUO, TWOEPO, C, CMl, PI, TOPI, PI TO, OS, OTW,  OTF 
COMMON  /PARFUN/  PHI , RHOQ, RHOQI , TRHOQ,  TRHOQI , RHOSQ,  RHOSQI , 03 
EQUIVALENCE  (RB3, BT3), (RB4, BT4) 

DATA  PC, PC2, PC3, ZR0/4*0.  /, EPS/1.  E-10/ 

C COMPUTE  DELTA-THETA 

03=1.  /3. 

PIM1=1. /PI 

TOP  I Ml  = 1.  /TOPI 

F0PIM1  = . 54LrOPIMl 

RCDTM1=RAD-»^CDTM1 

TPCTM1=RCDTM1*T0PIM1 

FPCTMl=RCDTMm-FOPIMl 

DTH=PI/(N1-1 ) 

SDTHM1  = 1.  /SIN(DTH) 

DTH2=DTH*.  5 

DTHM1  = 1.  /DTH 

DTHM1H=F0PIM1/DTH 

DTHM1C=DTHM1H*RCDTM1 

SDTH2=SIN(DTH2) 

SDT2M1  = 1.  /SDTH2 

C COMPUTE  NUMBER  OF  PATCHES  IN  STRIP 

TH=PI 
IS=2 

NPHI ( 1 )=1 
NPHI (Nl )=1 
N1H=(N1+1 ) /2 
DPH=1. 

IF( IPATCH.  EQ.  0)  THEN 

C NUMBER  OF  PATCHES  IN  A STRIP  HAS  TO  BE  A POWER  OF  2 

NPHX=1 

90  NPHX=NPHX*2 

1F(NPHX.  LT.  N3)  GO  TO  90 
DO  100  1=2. NIH 
TH=TH-DTH 
STH=SIN( TH) 

NPH=(N2-N3)^^STH+N3+.  5 
IF(NPH.  GT.  NPHX)  NPHX=NPHX<^2 
NPHI ( I )=NPHX 
IS=IS+NPHX 

DPH=AMIN1 (DPH,  STH*SIN(PI/NPHX) ) 

100  CONTINUE 

ELSE 

DO  105  1=2, NIH 
TH=TH-DTH 
STH=STN(TH) 
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NPH=(N2-N3)*BTH+N3+.  5 
NPHI ( I )=NPH 
IS=IS+NPH 

DPH=AMIN1 <DPH,  STH#SIN(PI/NPH) ) 

CONTINUE 
END  IF 

TO  ASSURE  SYMMETRY 
IF(M0D(N1,  2).  EQ.  1 ) N1H=N1H-1 
DO  1 10  J=I, Nl-1 
NPH=NPHI <N1H) 

NPHI < J)=NPH 
IS-=IS+NPH 
N1H=N1H-1 
CONTINUE 

CHOOSE  AN  ADEQUATE  (MINIMUM)  TIME  INTERVAL 
DT1  = 1.  6*RAD*CM1*AMIN1 (SDTH2,  DPH) 

IF(DT.  EQ.  0.  ) THEN 

USE  DTI  IF  NO  TIME  INCREMENT  IS  GIVEN 
DT^^DTl 
ELSE 

IF(DT1.  LT.  DT)  THEN 

USE  DTI  IF  GIVEN  TIME  INCREMENT  IS  TOO  LARGE 
PRINT  3,  DT, DTI 
NI=NI*DT/DT1 
DT=DT1 
END  IF 
END  IF 

PRINT  THE  COMPUTED  VALUES 

FILE  'OUTPT'  IS  USED  BECAUSE  THE  OUTPUT  FILE  TENDS  TO  FILL  UP 
IF( IDEB.  EQ.  1 ) THEN 

LENGTH- I S* 1 33/ 1 024+ 1 0 

CALL  Q5RQUEST(  'LFN=^  'OUTPT',  'MXR=',  133,  'LEN=^ LENGTH ) 

OPEN <7, FILE= 'OUTPT ' ) 

WRITE(7,  ' ( IHl ) ' ) 

WRITE(7,-}f)  'NUMBER  OF  PATCHES  IN  STRIPS' 

DO  120  1=^1,  Nl,  10 

WRITE (7.  1 ) NPHI ( I ) , NPHI < I + l ),  NPHI (1+2) , NPHI ( 1+3) , NPHI ( 1+4) , 

1 NPHI  (1+5),  NPHI  (1+6),  NPHI  (1+7),  NPHI  ( I +-8) , NPHI  (1+9) 

CONTINUE 
WRITE (7,  '(//)') 

END  IF 

PRINT  'NUMBER  OF  PATCHES  IS  ',IS,  ',  NUMBER  OF  ROWS  IS  ',LL 
IF( IS.  GT.  LL)  STOP  11 
IS1=^IS-NPH 
rn=p  I 

Bl=2.  *(RAD«SIN(DTH*.  25)  )^t*2 
IN=1 

COMPUTE  COMPONENT  OF  RADIUS  VECTOR,  SURFACE  ELEMENT,  TANGENTIAL  UNIT 
VECTORS,  ETC. 

BOTTOM  CAP,  THETA  ==  PI 
XRl ( 1 )=0. 

XR2( 1 )=0. 

XR3( 1 )=-RAD 
SRF(  1 )==S1 
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TF( 1 )=0. 

THETA  ( 1 ) =^P  I 
PHI  I ( 1 )~0. 

DPMI  I ( .1  ) -TOPI 
THl ( 1 )"1. 

TH2(  1 )"0. 

TH3(  1 ):=0. 

PHI  < .1  )=0. 

PH2(1 )=1. 

NP2=.l 

NP3-^NPHI  (2) 

I MB 2=1 

INB3=2 

DPHI2--=T0PI 

DPHI3=TDPI/NP3 

IF( ICORR.  NE.  0)  THEN 

C COMPUTE  COEFFICIENTS  FDR  SELF-PATCH  CORRECTIONS  FOR  CAPS 

C0<  1 ) = 1,  +.  25^5-DTH 
D0(  1 )=.  062!5*DTH-»-*2*RCDTMl 
D0(  IS)-=.  03125^<-DTH^^RCDTM1 
CALL  UD(PI. DPHI3, NP3,  1. 2, NL, NR, FCl,  FC2) 

CALL  UD(PITO, DPMI 3, NP3,  1, 2. NUl, NU2,  FUJ  , FU2) 

CALL  UD<PIT0«3.  , DPHI3, NP3,  1 , 2, ND 1 , ND2,  FD 1 , FD2 ) 

END  IF 

C REGULAR  PATCHES 

DO  200  1=2, Ml -1 
TH=TH-DTH 
STH=SIN(TH) 

CTH=COS(TH) 

'RSTH=RAD^P3TH 
RCTH=RAD-»tCTH 
T F I =CM  ( R AD+RCTH  ). 

NP1=NP2 
NP2=NP3 
NP3=NPHI  ( I-f  1 ) 

INB1=INB2 

INB2=INB3 

INB3=INB2+NP2 

DPHI 1=DPHI2 

DPHI2=DPHI3 

DPHI3=T0PI/NP3 

PHI=0. 

C 

C SOLID  ANGLE  FOR  THE  PATCH 

C 

SOLG=P  I M 1 ^^DPH 1 2#STH*SDTH2*R  AD^^^^2 

C 

C PATCH  COEFFICIENT 

C 

IF (ICORR.  NF.  0)  THEN 

C COMPUTE  COEFFICIENTS  FOR  SELF-PATCH  CORRECTIONS 

STH,DPH=STH-«-DPH  1 2 
STHDPH1=1. /STHDPH 
DGAMA2=DTH**?+STHDPH**2 
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DGAMMA=:SGRT  ( DGAMA2 ) 

AL 1 =ALOG  < ( DGAMMA+DTH ) / ( DGAMMA-DT  H ) ) 

AL2=AL0G< (DGAMMA+STHDPH) / (DGAMMA-STHDPH) ) 

PC2=F0P  I M 1 -»-STHDPH*AL  1 
PC 3=F0P  I M i }«-DTH-»-AL2 
PC-^PC2-PC3 

ATDG=,  5*DGAMA2*ATAN<STHDPH*DTHMi ) 

TAB=.  5*DTH*STHDPH 
BNPH2=.  5»fPIT0*STHDPH^^^^-2 
PC7=FPCTM  1 < TAB+ATDG-3NPH2  ) 

PCS=FPCTM  1 -»■  < TAB-ATDG^SNPH2  ) 

PD=PC7-PCS 
END  IF 

DO  190  J=i,NP2 
I 1 N+ 1 
CPH=--COB(PHI  ) 

SPH^SIN<PHI ) 

XRU  IN)=RSTH-^CPH 
XR2(  IN.P^RSTHr^SPH 
XR3(  IN)^=RCTH 
SRF( IN)=SOLG 
THETA ( IN)=TH 
PHI  I ( IN)=PHI 
DPHI I ( IM)  -=^DPHI2 
THl  ( IN):=CTH->^^CPH 
TH2<  IN)=CTH^^-BPH 
TH3-;  IN)==-STH 
PHI  ( IN)=-^-SPH 
PH2  (IN)  ==CPH 
TF( IN)=TFI 
IF( ICGRR.  NE.  0)  THEN 
CO(IN)=^PC 
D0( IN)=PD 
IF(J.  EG.  1)  THEN 
NL( IN)=IN+NP2-1 
ELSE 

NL  ( I N ) :=  I N- 1 
END  IF 

IF(J.  EG.  NP2)  THEN 
NR ( IN)=IN-NP2+1 
ELSE 

NR ( IN)^IN+1 
END  IF 

CALL  UD(PHI,  DPHIl,  NP  1 , IN,  INBl,  NDl,.  ND2,  FDl,  FD2) 
CALL  UD(PHI, DPHI3, NP3,  IN,  INB3,  NUl, NU2,  FUl, FU2) 
END  IF 

PHI=PHI+DPHI2 
190  CONTINUE 
200  CONTINUE 
C TOP  CAP,  THETA  = 0 

XRl  < IS)=^0. 

XR2( IS)=0. 

XR3( IS)=RAD 
SP.F(  IS):^S1 


63 


THETA(  IS):^0. 

PHI  I ( IS)  =0. 

DPHU  ( IS)=TOP:t 
THl ( IS)"1. 

TH2(  IS)=--^0. 

TH3(  IS)==--0, 

IF< ICORR.  NE,  0)  THEN 

CALL  UD  (P.T,  DPHI2,  NP2,  IS,  IND2,  NL,  NR,  PHI,  PH2) 

FCl (2)=PH1 ( IS) 

FC2(2)=PH2( IS) 

CALL  UD(PI  TO,  DPHI2,  NP2,  IS,  INB2, NUl, NU2,  FUl,  FU2 ) 
CALL  UD(PITG*3,  DPHI2,  NP2,  IS,  INB2, NDl, ND2, FDl, FD2 ) 
END  IF 
PHI  ( IS)=-0. 

PH2(  ,[S)^='-1. 

TF(  IS)^::-2.  i^CM,1,*RAD 

C CONFUTE  QUANTITIES  RELATED  TO  TWO  PATCHES 

DLl  = DTH-«-RAD 
NMAX^=-0 
I N==  1 
NPNT^-“-0 
NP:-..=0 

DO  240  1=^1,  IS 
X-^^^XRl  ( I ) 

Y=--XR2  ( I ) 

Z=--^KR3(  I ) 

IH~THETA( I ) 

PHI=PHI  Id) 

STi-L^SINdH) 

RSTH=RAD-J>STH 

NM=1 

IF(NP.  EQ.  0)  THEN 
NPH^^^NPHI  (IN) 

DPH^TDPI/NPH 

DL2--=DPH*RSTH 

rMAX=SQRT(2  )-i^AMAXl  (DLl,  DL2) 

IFd.  E(T  1.  OR.  I.  EQ.  IS)  RNAX=,1.  5->^DTH<^RAD 
NP--^NPH 
I N--'  I N+ 1 
END  IF 
NP--=NP-t 

DRl  (1;  1S)  = X-XRI  ( 1.;  IS) 

DR2(  ii  IS)  =-Y-XR2(  li  IS) 

DR3di  IS)^Z-XR3(  1;  IS) 

DR(  1;  IS)“V3QRT(DR1  ( li  I S ) -i^-t'-;2+DR2 ( li  I S ) -«-«-2-^DR3 ( 1;  IS) 
1 DR  d ; I S ) ) 

DRN3(  1;  JS)=CDTNH^DR  di  IS) 

IND( li IS)=VINT (DRN3( 1; IS)i IND( 1, IS) ) 

C DETERMINE  WHERE  SHIFTS  IN  ITERPOLATION  ARE  NEEDED 

BT3(  1,  I;  IS)=^INDdi  IS).  EQ.  2 
R T 4 d , I , I S ) I ND  d i IS).  EQ.  1 
IFdCORR,  LT.  2)  GO  TO  230 

C COMPUTE  COEFFICIENTS  FOR  NEIGHBORING-PATCH  INTEGRALS 

DO  220  J=l, IS 
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TF(  J.  £G.  I.  DR.  DR(  J».  GT.  RMAX)  GO  TO  220 

IF(NM  GT.  LL20  OR . NPNT.  GT.  LL 1 2 ) THFN 

PRINT  NPNT-- ' , NPNT,  ' , J-=  ' , J 

STOP  12 
END  IF 
NE I < I , NM ) 

NPNT=NPNT-H 
NPOINT(  I,  Nh)==NPNT 
IF(J.  EQ.  1.  OR.  J.  EQ.  IS)  THEN 
NEIGHBORING  PATCH  IS  A SPHERICAL  CAP 
RH0Q=^STH^SDT2M1 
TRH0Q^2  *RHDQ 
RHnQH=STH>5F0PIMl*SDTHMl 
RH0QCT=RH0GH^^RCDTM1 
R HOG!  1 = 1.  /RHOQ 
TRH0QI=2.  frRHOOI 
RH0SQ=1.  f-RH0Q^t*2 
RH0SQI  = 1.  +RH0G.I-S-S-2 
RHCTm=STH■«-RCDTNi^^TOPIMl 
RHCTM1H=.  5*RHCTN1 


USE  QIDB  FROM  CMLIB  TO  DO  NUMERICAL  INTEGRATIONS 

CP003  INCLUDES  AN  ADDITIONAL  FACTOR  RHQO 

CALL  0tDB<F003, ZRO, TOPI, EPS, CP003, ER, KF, I FLAG) 

IF(  IFLAG.  GT.  3)  THEN 

PRINT  'PROBLEM  IN  CP003=  ' , CP003,  : ERROR=  ',ER, 

1 NO.  OF  COMPUTATIONS^ KF,  ',  FL.AG=',  IFLAG 

STOP 
END  IF 

CALL  QIDB (F103,  ZRO,  TOPI, EPS, CP  103,  ER, KF,  IFLAG) 

IF (IFLAG,  GT.  3)  THEN 

PRINT  'PROBLEM  IN  CP  1 03=  ' , CP  1 03,  ' ERROR=',ER, 
1 ',  NO  OF  COMPUTATIONS^',  KF,  ',  FLAG=  ',  IFLAG 

STOP 
END  IF 

CALL  Q1DE(F013,  ZRO..  TOPI,  EPS,  CP013.  ER,  KF,  IFLAG; 

IF( IFLAG  GT.  3)  THEN 

PRINT  'PROBLEM  IN  CP013= ',  CP013,  ',  ERROR=  ',ER, 
1 ',  NO.  OF  COMPUTATIONS^ ',  KF,  FL  AG= ',  IFLAG 

STOP 
END  IF 

CP 203  INCLUDES  AN  ADDITIONAL  FACTOR  1/RHOO 

CALL  Q1DB(F203.  ZRO, TOPI, EPS, CP203, ER, KF, IFLAG » 

IF (IFLAG.  GT.  3)  THEN 

PRINT  'PROBLEM  IN  CP203=  ' , CP203,  ' , ERROR=',ER, 
1 ',  NO  OF  COMPUTATIONS-^',  KF,  ',  FL.AG= ',  IFLAG 

STOP 
END  IF 

CPI  13  INCLUDES  AN  ADDITIONAL  FACTOR  1/RHOO 

CALL  QIDB (FI  13,  ZRO, TOPI, EPS, CPI  13, ER,  KF,  IFLAG) 

IF( IFL AG.  GT.  3)  THEN 

PRINT  'PROBLEM  IN  CP  11  3=  ' , CP  1 1 3,  ',  ERROR=',ER. 
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1 NO  OF  COMPUTATIONS=--=S  KF,  S FLAG--=  ^ IFLAG 

STOP 
END  IF 

CP023  INCLUDF:,S  an  additional  FACTOR  1/RHOO 

CALL  QIDB (F023,  ZRO, TOPI, EPS,  CP023,  ER,  KF.  IFLAG) 
IFdFLAG.  GT.  3)  THEN 

PRINT  'PROBLEN  IN  CP023«  ' , CP023,  ' i ERROR=',ER, 
.1  L.  NO  OF  CONPUTATIONS--=- KF,  FLAG= I FLAG 

STOP 
END  IF 

CALL  GtDB<F002, ZRO, TOPI, EPS, CP002, ER, KF , IFLAG) 
IFdFLAG.  QT.  3)  THEN 

PRINT  ■«-,  'PROBLEM  IN  CP002=-' , CP002,  ' ; ERROR^',ER, 
I NO  OF  COMPUTATIONS= ',  KF,  FLAG-',  IFLAG 

Slop 
END  IF 

CP. 1.02  INCLUDES  AN  ADDITIONAL  FACTOR  1/RHOO 

CALL  Q.1DB-:F1G2,  ZRO,  TOPI,  EPS,  CP102,  ER,  KF,  IFLAG) 

IF(  IFLAG.  GT.  3)  THEN 

PRINT  'PROBLEM  IN  CP  .1 02=  ' , CP  .1  02,  ' ; ERROR- ',  F:R  , 

i ',  NO  OF  COMPUTATIONS- ',  KF,  ',  FLAG-',  IFLAG 

STOP 
END  IF 

CP012  INCLUDES  AN  ADDITIONAL  FACTOR  1/RHOO 

CALL  Q1DB(F012,  ZRO,  TOPI, EPS, CP012, ER, KF,  IFLAG) 

IF< IFLAG  GT.  3)  THEN 

PRINT  'PROBLEM  IN  CPO 12=  ' , CPOl  2,  ' , ERROR- ',ER, 
',  NO  OF  COMPUTATIONS- KF,  FLAG-',  IFLAG 
STOP 
END  IF 

CP202  INCLUDES  AN  ADDITIONAL  FACTOR  l/RH00^<-«2 

CALL  G;lDB(F202,  ZRO,  TOPI,  EPS,  CP202,  ER,  KF.  IFLAG) 
IFdFLAG.  GT.  3)  THEN 

PRINT  -ti-,  'PROBLEM  IN  CP202- ' , CP202,  ' ; ERROR-', ER, 
1 NO  OF  COMPUTATIONS- ',  KF,  ',  FLAG=',IFIAC 

STOP 
END  IF 

CPI  12  INCLUDES  AN  ADDITIONAL  FACTOR  1/RH004»h^2 

CALL  G1DB(F112, ZRO, TOPI, EPS, CPJ 12, ER, KF, IFLAG) 
IFdFLAG.  GT.  3)  THEN 

PRINT  'PROBLEM  IN  CP  1 1 2- ' , CP  d 2,  ' , ERROR-', ER, 
..  NO  OF  COMPUTATIONS- ',  KF,  ',  FLAG- ',  I FLAG 
STOP 
END  IF 

CP022  INCLUDES  AN  ADDITIONAL  FACTOR  I/RHOOm-4^2 

CALL  GIDB (F022, ZRO, TOPI, EPS, CP022. ER, KF, IFLAG) 
IFdFLAG  GT.  3)  THEN 

PRINT  'PROBLEM  IN  CP022- ' , CP022,  ' , ERROR- ',ER, 
i ',  NO  OF  COMPUTATIONS- ',  KF,  ',  FLAG-'.  IFLAG 

STOP 
END  IF 
CPH=COS<PHI ) 

SPH-SIN(PHI ) 

CO  J ( 1 , NPNT  ) -TOP  I M 1 ■«•  ( CPH-K-CP003-  CP  1 03  ) 
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COJ  ( 2.  NPNT  ) =TOP  I Ml  * < SPH-J^CP003-CP0 1 3 ) 

COJT( 1, NPNT )=RHCTMH^<CPH*CP002-CP 102) 

CO JT ( 2. NPNT ) =RHCTMi * ( SPH*CP002-CP01 2 ) 

IF(  J.  EQ.  1 ) THEN 

COJ  ( 3 . NPN  r ) =FOP  I M 1 -K-S 7 H*  ( CP003-CP203--CP023 ) 
C0JT(3;  NPNT)=RHCTM1H-J^STH*<CP002-CP202-CP022> 
ELSE 

COJ  ( 3,  NPNT  ) =-FOP  I Ml-K-STH*  ( CP003-CP2C3-CP023  ) 

CO JT  ( 3,  NPNT  ) = -RHCTM 1 H^^STH-*^  ( CP002-CP202-CP022 ) 
END  IF 

COJU(  1,  NPNT)  = RH0QH-k-<CPH-ij-CP103-CP203) 

C0JU(2,  NPNT)=RH0QH^^(SPH*CP103-CP1  13) 

COJVd.  NPNT)=RH0QH*<CPH*CP013-CP1  13) 

COJ  J < 2,  NPNT  ) =^RHOQH-{?-  < SPH-»^CPO  1 3-CP023  ) 

COJTU<  1,  NPNT)=RH0QCT^^(CPH#CF102-CP202) 

COJTU  ( 2,  NPNT  ) =^^RHOQCT-f^  ( SPH-«-CP  102--CP  112) 

COJTV(  1.  NPNT)=^RH0QCT-»-(CPH#CP012--CPl  12) 

CO  JT  J (2.  NPNT  ) =RHOQCT*  < SPH-«-CP01 2-CP022  ) 

ELSE 

C OTHER  NEIGHBORING  PATCHES 

THP=--THETA(  J) 

PHIF--PHI I ( J) 

GTHP^^^SINCI  HP  ) 

CTHP*---COS  ( THP  ) 

CSCTHP=^1,  /STHP 
C0TTHP=CTHP*C3CTHP 
BPHP---SIN(PHIP  ) 

CPi  !P-=COS(PHIP  ) 

STCP==STHP-«-CPHP 
STSP- ST HP ^SP HP 
CTCP=CTHPii-CPHP 
CTSP=--CTHPi'cSPHP 
DELTH=THF-TH 
DELPHI=PHIP-PHI 
IF(I.  EQ.  1.  OR.  I.  EQ.  13)  THEN 
DEEP H 1=0. 

ELSE. 

IF(DEi-PHI.  GT  PI  ) DEl..PHI=DELPHl-TOPI 
1F(DE!.PHT  LT.  -PI)  DELPh  I =DELPH  I +TOP  I 
END  IF 

STD£LP=STHP-)tDELPHI 
THTP=DEI..TH+.  5*DTH 
THTP2=THTP-«->^2 
TH  TP3=THTP2-«-THTP 
THTP4=THTP3-«-THTP 
THTM=DEI.TH-.  5*DTH 
THTM2-=THTM«-»-2 
THTM3=THTM2*THTM 
THTM4=THTM3*THTM 
DPHI=DPHI I ( J) 

PHTP=(  DELPHI +.  S^^DPHI  ) i^STHP 
PHTP2=PHTP^(•■«■2 
PHTP3=PHTP2-R-PHTP 
PHTP4-PHrP3*PHTP 
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PHTM=M DELPHI-,  S^^DPHI  > *STHP 

PHTM2-PHThH<-*2 

PHTh3=PHTM2*PHTM 

PHTM4-PHTM3*PHTM 

GAM 1 T~SGRT ( THTM2+PHTM2 ) 

GAM1TI  = 1.  /QAMIT 
GAMlTI2--=GAMlT.T*->^2 
GAM  1 T 1 3-GAM 1 T 1 2^«•GAM  1 T I 
GAM2T^^SQRT  ( THTP2t-PHTM2 ) 

GAM2Ti  ==l.  /GAM2T 
GaN2T  T 2=GAM2T  I >t^2 
GAM2T .[  3---GAM2  V 1 2^^GAM2T  I 
GAM3T  ^^SGRT  ( THTP24-R|f[ P2  ) 

GAM3T:i^  i.  /GAM3T 
QAM3T  I 2=^-GAM3T .( **2 
GAM3T  ;i  3^G AM3T  1 2-t^GAM3T  I 
GAM4T-:SGRT  ( THTM2+PH  f P2  ) 

GAM4T1--1.  /GAM4T 

GAM4  f T2^GAM4TI-»-h2 

GAM4T  :i  3 =GAM4T  1 2-R-GAM4  T I 

THMAX2==AMAX1  (THTM2.  THTP2)*1.  E-J2 

PHMAX2--^AMAX1  (PH7M2,  PHTP2.>-k-1.  E-12 

IF(THTM2.  LT.  PHMAX2.  AND.  PHTP.  LT.  O ) THEN 

C103=--AL0G(PHTM-«  (GAM2T+PHTM)/<  ( GAM3T-^PHTP  ) »-PHTP)  ) 
ELSE  IF(THTP2.  LT.  PHMAX2.  AND.  PHTP.  LT.  0.  ) THEN 

C103--=  AL0G(  (GAM4T+PHTP)i<-PHTP/(PHTMi^(GANl  T-t  PHTM/  ) ) 

El...  BP 

C 1 03^AL0G  ( ( GAM4T  + PHTP  ) * ( GAM2T  * PHTM  ) / ( v GArl3T-4  PHTP  ) fr 

t (GAhiT^-PHTM)  ) ) 

END  IF 

IF(PHTM2.  LT  THMAX2.  AND.  THTP.  LT.  O.  ) THEN 

C013=-AL0G(THTM-k-(GAM4T<-THTM)/  < (GAN3T+THTP  )^<THTP  > > 
ELSE  IF  (PHTP2.  LT.  THMAX2.  AND.  THTP.  LT  0.  ) THEM 

C013'^ALGG(  (GAM2T4-THTP  ) b-THTP  / ( THTM*  ( GAM  1 T + THTM  ) > ) 
ELSE 

C013=^AL0G(  <GAM2T^-THTP)*(GAM4T  +-THTM)/(  <GAN3TfTHTP)  >* 

•i  (GANi  T-t-THTM)  ,>  ) 

END  IF 

IF(PHTM,2.  LT.  THMAX2)  THEN 

C203-=  -PHTM*ABS  ( ALOG  ( THTP  /THTM  ) ) 

2 5-sPHTP*AL0G<  (GAM3T+THTP  )*(GAM4T-THTM)  / 

3 ( ( GAN3T-Tt-n  P ) * ( GAM4T  +THTM ) ) ) 

ELSE  IF(PHTP2.  LT.  THMAX2)  THEN 

C203=-.  5*PHTM-«-ALDG(  ( GAN2T -t-^THTP  )»>■(  GAN  1 T -THTM  ) / 

1 ( ( GAN27  -THTP  ) < GAM  1 T+THT  M ) ) ) 

2 ^-PHTP^^■ABS<  ALOG  (THTP /THTM)  ) 

ELSE 

C203=-^-.  D-«-PHTM-frALOG<  < GAM2T  + THTP  ) * ( GAM  1 T-TH  I M ) / 

1 ( (GAN2T-THTP )*(GAM1T+THTM) ) ) 

2 +.  5^PHTP-»^AL0G(  <GAM3T-(-THTP)-»KGAM4T-THTM)/ 

3 ( ( GAN3T -THTP ) * ( GAM4T+THTM ) ) ) 

END  IF 

C 1 1 3=-GAMl T+GAM2T-GAM3T+GAM4T 
IF(THTM2.  I T.  PHMAX2)  THEN 
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C023= - T HTM*ABS ( ALGO ( PHTP/PHTM  > ) 

2 4 . 5*THTP*ALGG  ( ( GAM3T4-PHTP  > * ( GAM2T-PHTM  ) / 

3 ( ( GAM3T- PHTP ) * < GAM2T4  PHTM ) ) ) 

ELSE  IF(THTP2.  LT.  PHMAX2)  THEN 

C023^-.  5*THTN*ALOG<  ( GAM4T +PHTP  ) ( GAMIT-PHTM ) / 

1 ( (GAN4T-PHTP )*(GAM1T+PHTM) ) ) 

2 4- THTP*ABS(ALOG(PHTP/PHTN)  ) 

ELSE 

C023^~  . S-t^THTN^ALOG  ( ( GAM4T4-PHTP  ) ( GAMIT-PHTM  ) / 

1 ( (GAi'-14T-PHTP  ) *(GANlT4-PHTIi)  ) ) 

2 4 . 54^THTP#AL  OG  ( ( GAM3T4  PHTP  ) # < GAM2T-PHTM  ) / 

3 ( ( GAM3T-  PHTP  ) ( GAM2  T + PHTM  ■>  ) ) 

END  ir 

C125=03-«  ( (PHTP4^(GAM3TI GAM4TI  ) -PhTM-s- ( GAM2T  I -GAM  i T I ) > 

1 4-C103) 

C035-^03k->.  (THTM4J-(GAM4TI  -GANITI  ) -TH  f P-!?- ( GAM3T I -GAM2T I ) 
4-2.  *COi3) 

C225:---^Q3-s-'.  fHTP*(PHTP*GAM3T  I -PHTPHJ-GAN2T I ) 

1 --(THTMiv(FHTP*GAI'14  ri-PHTri*GAMlT  1 ) ) ) 

C 1 35  = THTM2-«'  ( 03-»-THTM2  r-  ( GAMl  T I 3-GAM4T 1 3 ) -GAN  1 T I 

1 4-GAN4TI  )-*-PHTP4*03<^(GAN4TI3-GAN3T.T3) 

2 (-THTP2  «■  ( 03*THTP2*  ( GAN3T  I 3-GAN2TI  3 ) -GAN37  I 

3 4-GAN2TI  )4PHTM4-fr03«-<GAM2TI3-GANlTI3) 

DTHT=  fHTP -THIN 

DPHT:-^PHTP-PHTN 

PHT  1 =^ATAN2  ( PHTM,  THIN  ) 

PHT'2:^-ATAN2  ( PHTN,  THTP  > 

PH  f3=ATAN2  ( PHTP THTP  ) 

PHT4=^ATAN2(PHTP,  THTN ) 

PHT.12==PHT1-PHT2 

PHT23==PHT2-PHT3 

PHT34=PHT3-PHT4 

PHT41=PHT4-PHT1 

IF(PHTM  LT.  0.  . AND.  PHTP.  GT.  0.  . AND.  THTP.  LT.  O.  ) THEN 
PHT23  =PHT234-TGP  I 
PHT4  i:=PHT41-T0PI 
END  IF 

ALG  .1  2=-ALGG  ( GAM  1 T>GAN2T  I ) 

ALG23=^AL0G  ( GAN2TR-GAN3T I ) 

ALG34=ALGG(GAN3T*GAN4TI ) 

ALG41=AL0G(GAN4T-»GAN1TI  ) 

C102=“THTPsPHT234-PHTP4S<^LG34-THTM*PHT4  1+PHTN*ALG12 
CO  1 2=^-PHTP  )*PHT34-TH1  P*ALG23-- PHTM*PHT  1 2-THTN#ALG4 1 
C202-=.  5-k-(DTHT-«-DPHT4PHTN2*PHT12-THTP2k^PHT23 
1 4PHTP2+:PHT34-THTM24iPHf41  ) 

C.l  12-  5t^(PHTM2-t^ALG12-THTP2-KALG23 
1 4 PHT  P2-K  ALG34“THTN24i  ALG4  1 ) 

€022=--.  5 ■:  DTHT*DPHT-PHTN2*PHT  1 2 < THTP2+^PHT23 
1 -PHTP2-i<PHT344-THTM2*PHT41  ) 

C 1 24  = . ( PHTN3*  < GAM2T 1 2-GAN 1 T 1 2 ) -THTP->s-PHT23 

1 HHTP24«  (PHTN*GAN2^I2-PHTP^^GAM3T  12) 

2 4 PHTP3  K-  ( GAM4T  1 2-GAN3T 1 2 ) -THTM*PHT4  1 

3 h-THTN2«  (PHTP*GAM4TI2-PHTM-«-GAMl  r 121  ) 

C034--.  PHTM*PHT12-PHTN2*(THTP*GAN2TI2-THTNkGAN1TI2) 
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1 >THTP-«-vPHTM2#0AM2TI2"PHTP2*GAM3TI2-2,  -»ALG23) 

2 - PHTP*PHT34-PHTP2* ( THTM^GAMATK’-THTP^GAMST 1 2 ) 

3 MHTI''|>^(PHTP2*GAN4TI2-  PHTM2*GAM1TI2~2.  •»-AL.G4l  ) ) 

C224^.  25*  ( ~PHTM2*PHT  1 2 + PHTM3*  ( THTP*GAM2T 1 2 --THTM*GAM  1 T 1 2 ) 

1 -THTP2*PHT23+THTP3*  < PHTM*GAH2T 1 2-PHTP*GAM3T 1 2 ) 

2 -PHTP2*PHT34+PHTP3*( rHTM*GAM^TI2-7HTP*GAM3TI2) 

3 --THTM2*PHT4 1 +THT  M3*  ( PHTP *GAM4T .1 2--PHTM*GAM  I T 1 2 ) ) 

Cl 34--=.  25*(PHTM4*(GAM2TI2-GAM1T  12) 

.1  i^THTP2*(PHTM2*GAM2TI2-"PHTP2*GAM3TT2-2.  *ALG23  ) 

2 +PHTP4*  ( GAM4T I 2-GAlM;3T  1 2 ) 

3 *THTM2*(PHTP2*GAM4TI2-PHTM2*GAMlT12  -2.  *ALG'-'V1  ) ) 

DTHCOT=DELTH*COTTHP 

CljETH-ClOS^i  < -1.  +DTHC0T)  + (-C203<-l  5 «»C225  ) *COTTHP 
1 -3.  *C125*DTHCOT 

C0£PH==C013*(  -1.  +2.  *DTHC0T>+CTHP«-DELPH1*C103 
1 *C0TTHP*(“2. *C113vl  5*C135)-3.  *DTHC0T*CG35 

COEIMO=.  5*C203--DE.LTH*C  1 03 
C 0 El  R . 5 »■  C 0 2 3 * C S C T H P - D E L P H I * C 0 1 3 
CTHTH=DELTH*C 1 03-C203 
CTHPH==dELTH*CO  1 3 -C  .1  1 3 
CPffrH-3TDELP*C  1 03-C  1 1 3 
C P HP  H^C>  T DEEP  *C  0 1 3 -C  02  3 

COJ  1:  JTHTH.  2; JTHPH  3 J1 HM.  4-JPHTH,  5 JPHPH  6 JPHN 

C O J < 1 .■  NP MT  ) = “TOP  I M 1 *C  1 PIP H 

COJ  ( 2,  NPNT  ) =TOP  IMl  * < COENO+COEF?*S T HP  + CTHTH  ) 

COJ<  3>  NPNT)=--T0P  .TM1  ( -CDEPH-t-CPHT  H*COT  T HP  ) 

C0J(4,  NPNT  >=T0PIMl-*^r(-C0EN0-C0FR>J-5THP-  CPHPH) 

CGJ(  5.  iMPNT)=TOPIMl*CPHTH 

CO  J <6..  MPNT  >=^T0PIMl*(C0ElH-t-C0ER*CTHP4-CPHPH*CnTTHn  ■ 

COJU'  .1,  MPNT)=-DTHM1H*C1  HPH 
CGJU<2, NPNT )=DTHM1H*CTHTH 
STDMIH^EDPIMI  / (STHPi^DPtHI  ) 

CGJV(  1..  NPNT)==-STDM1H*CPHPH 
CO JV ( 2>  NPNT ) =STDM1 H*CPHTH 

C0ETH=C102*( -1  +DTHC0T)  + ( -C202«-C22^  ?*CCn  I HP 
.1  --2,  *C124*DTHC0T 

C0EPH==C012*( -1  -t-2  *DrHCOT  )+CTHP*DELPHl  + C 102 
] -i-CGl  T HP*  < -2.  *C1  l2*C134)-2.  « DTHCDl  hC034 

CQt-.NO=-^.  5*C202--DELTH*C102 
COER=.  5*C022*CSCTHP -DELPHI *C012 
CTHTH=DELTH*C 102-C202 
CTHPH=DELTH*C012-C1 12 
CPHTH=^STDELP*C102-C1  12 
CPHPH=S  rDELP*CO  1 2--C022 
COJTd,  NPNT)=-TPCTM1*CTHPH 

CO  JT  ( 2,  NPNT  ) =TPCTMi  * ( COENO+COER  *STHP-t  CTH  fH  ) 

CO  JT  < 3,  NPNT  ) =TPCTM1  ft  < -COEPH-t-CPHTH*COT  IMP  ) 

CU  JT  (4..  NPNT  ) =TPCTM1  * ( -COENO-COER*STHP-CPHPH  ) 

COJT( 5. NPNT)=TPCTM1*CPHTH 

C0JT(6,  NPNT  )=^^TPCTMl  ft(COETH+COERftCTHP>CPHPHftCOJTHP  > 

COJTU( 1 , NPNT)=-DTHM1C*CTHPH 
CGJTU(2/  NPNT)=-DTHIilC*CTHTH 
STDM1C=STDM1H*CRDTM1 
COJTV( I,  NPNT)=-3TDM1C*CPHPH 
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C0JTV(2, NPNT)=STDM1C#CPHTH 
END  IF 
2S0  CONTINUE 

230  NEI(I,.  1)=NN 

I F ( NM . QT . NMA X ) NNA X = NM 
240  CONTINUE 

PRINT  ^NUMBER  OF  COLUMNS  NEEDED  IS  UNMAX-1,  ^ OUT  OF  19, 

1 ' TOTAL  NUMBER  OF  NEIGHBORING  PATCHES  IS  U NPNT, 

2 U.  W I TH  A hi  A X I MUM  OF  •- , LL 1 2 
RETURN 

1 FORMAT^  10(2X,  I 10)  ) 

3  FORMAT  GIVEN  DT-=UE9.  2,  ' GREATER  THAN  "MINIMUM"  DT=^UE9.  2) 

END 

SUBROUTINE  UB(PHI, DPHI, NP,  IN,  INB, NXl,  NX2,  FX1  , FX2) 

C 

C THIS  SUBROUTINE  IS  USED  TO  COMPUTE  QUANTITIES  NEEDED  FOR  PHI -INTERPOLATION 

C 

PARAMETER  (LL==L LG,  LC-LCG) 

D I MENS I ON  NX  1 ( LL ) , NX2  < LL ) , F X 1 < LL ) , FX2  < LL ) 

FI -PH I /DPHI 
Ml -FI 

Fi-Fl-REAL(N.l  ) 

NXl  ( IN)  = IMB-H41 
IFUMl  EQ  NP-  i ) THEN 
NX2‘:  IN  ^ -INB 
ELSE 

NX2( IN)=INB+N1+1 
END  IE 

IFUMi  . EQ.  N2)  THEN 
C FOFv'  POLAR  CAPS 

FXl ( IN) -COS (PHI ) 

FX2( IN )- SIN (PHI  ) 

ELSE 

FXl ( IN)-1. -FI 
FX2( IN) -FI 
END  IF 
RETT'RN 

END 

r 

FUNCTION  FOOCKX) 

common  /PARFUN/  PHI , RHOQ, RHOQI,  TRHOQ.  TRHOQI , RHOSQ,  RHOSQI , 03 
C03X-C0S( X ) 

EKS-SQRT  ■;  RHOSG-TRHOQ-frCOSX  ) 

F003- 1 . / <;  EKS-<5  ( EKB  vRHOQ-COSX  ) ) 

RETURN 

END 

r 

FUNCTION  F103(X) 

COMMON  /PARFUN/  PHI, RHDQ, RHOQI, TRHOQ,  TRHOQI,  RHOSQ,  RHOSQI,  03 
CQGX=CGS(X) 

RCQGX-RHOQB-COSX 
TCGSX-2.  -«-RCOSX 
EKS-SQRT ( RHOSG-TCOSX ) 
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ECOSX=^EKS#CC)GX 
IF<COSX,  GT.  . 5)  THEN 

F"103=-^C0S(X+PHI  )-5^(  ALOG(  < RHOQ+RCOSX  ) / < EKS-^RCOSX-1 . ) > 4- 

1 2.  *COSX/ (EKSXEKS+RHOQ-COSX)  )- 

2 (1  -TCOSX)/(EKS*( 1. -ECOSX-RCOSX) ) ) 

ELSE  IF(CGSX.  GT,  5)  THEN 

Fi03=^C0S<X4PHI  )*(ALOG'  ( RHGQ+RCOSX  ) / < EKS+RCOSX“l . ) ) + 

1 2.  -frCOSX/  (EKS^^(EKS4  RH0Q~C0SX  ) > -- 

2 (1.  ^-ECOSX--RCOSX)/(EKS^<•(  1.  -C0SX**2>  ) ) 

ELSE 

Fi03=-C0S(  X4PHI  )*<  ALOG(  ( EKS-RCOSX  + 1 . ) / <RH0Q-RC05X  ) ) + 

1 2 )^COSX/ (EK5*(EKS+RH0Q-CC)3X  ) )-■ 

2 (1.  -TCOSX  ) / (EKS^^(  i -ECDSX-RCOSX  ) ) ) 

END  IF 

RETURN 

END 

F'JNCTKDN  F0i3(X) 

COMMON  /PARFUN/  PHI, RHOQ, RHOQI, TRHOQ, TRHOQI, RHOSG. RHDSGI, 03 
COSX^COSC  X .) 

RC0SX=:RH0Q-»^C0'3X 
TCDSX=^2.  -»-RC.OSX 
£KS==SQRT(RHOSG-TCOBX  ) 

ECOSX--EKS*CGSX 
IF(COSX.  GT.  f>)  THEN 

FO  13=^3  IN  (X+-PHI  (ALOG<  ( RHOQ-^-RCOSX  ) / ( EKS+RCOSX  - 1 ) ) + 

.1.  2.  -iiCOSX  / ( EKS^^-  < EKS+RHOQ-COSX  ) ) - 

2 (J,  "TCOSX)/(EKS-(^(  1.  -ECOSX-RCOSX  ) ) ) 

ELSF  IFCCOSX.  GT  - 5)  THEN 

F0i3=^SIN(  X^•PHI  ':AL0G<  (RHOQ+RCOSX)/ (EKS4RC03X-I.  ) )4 

1 2,  -«^COSX/vEKS«-(EKS+RHOQ~COSX)  )- 

2 (.1.  4EC0SX-RC0BX)/(EKS*v1.  -C0SX*-«2))) 

ELBE: 

F013=SIN<  X+PHI  )*(  AL0G(  ( EKS-RCOSX  + 1 . ) / ( RHOQ-RCOSX  ) )■♦ 

1 2,  -fs-COSX/ (EKS*<EKS+RHOQ-COSX  ) )- 

2 (1  -TCOGX ) / ( EKS# ( 1 . -ECDSX-RCOSX ) ) ) 

END  IF 

RETURN 

END 

FUNCTION  F203(X) 

COMMON  /PARFUN/  PHI. RHOQ. RHOQI, TRHOQ. TRHOQI. RHOSQ, RHOSQI. 03 
COGX=COS( X ) 

RCOSX=RHOQ*CC)SX 
TC0SX=2.  ^RCOGX 
EKS^--SQRT  ( RH0GQ-TC03X  ) 

ECOSX=EKS*CGSX 
IF(COSX.  GT.  . 5)  THEN 

F203-C0S(X+PHI  )*-«2-«-(RH0QI*EKS-1.  + 

1 3.  -«-COSX*ALOG(  (RHOQ-t-RCOSX)/(EKS+RCOSX-l.  ) ) 4- 

2 (4.  •s-C0SX**2-l . )/ (EKS-k-(EKG4-RH0Q-C0SX)  )- 

3 2,  -«-COSX*(  1.  -TC0SX)/<EKS-J^(  1.  -ECOGX-RCOGX  ) ) ) 

ELBE  IF(COSX.  GT.  5)  THEN 

F203--'-^C0B(X4  PHI  ) *^<2*  ( RHOQI *EKS- 1 . 4- 
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1 3.  -itCOSX*ALOG(  (RHOQ-»-RCOSX)/(EKS+RCOSX-l.  ))  + 

2 <4.  *C0SX**2~1.  )/<EKS*<EKG+RHOQ--COSX)  )- 

3 2,  5tC0SX*  ( 1 . +ECOSX-RCOSX  ) / < EKS*  < 1 . “C0SX^t*2 ) ) ) 

ELSE 

F203=C0S(  X+PHI  )*-«2#(RH0QHtEKS-l.  + 

1 3.  ■H-CDSX*ALOG(  (EKS-RC03X  + 1.  ) / ( RHOQ-RCOSX  ) ) + 

2 (4.  ■»^C0SX**2--1.  ) / <EKS-?<-<EKS+RHOQ“COSX  ) )- 

3 2 *COBX*<  1.  -TCOSX)/(EKS<^(  1.  -ECOSX-RCOSX  ) ) ) 

END  IF 

RETURN 

END 

FUNCTION  F;ll3(X) 

COMMON  /PARFUN/  PHI, RHOQ, RHOQI>  TRHOQ,  TRHOQI,  RHOSQ>  RHOSQI, 03 
CDSX--COS(  X ) 

RCOSX=RHOQ>tCOSX 
TCG5X=2.  *RCOSX 
EKS==SGRT>:RH0SQ  -TCOSX  ) 

EC05X=EKS^C0GX 
IF^COSX.  GT.  , 5)  THEN 

FI  13=C03(X-^PHI  )*5IN<X+PHI  )*(RH0QI-frEK3“l.  + 

1 3 -»^COSX*ALOG<  <RH0Q-i-RC0SX)/(EKS4-RC0SX-l.  ))  + 

2 (4.  ^COSX-s-K-2-1.  ) / (EKS*<EKG4-RH0Q“C0SX  ) )~ 

3 2.  -«C0SX*(1.  -TCOSX)/<EKS*a.  “EC0SX-RCQ3X ) ) ) 
else  IF(C0SX.  GT.  5)  THEN 

FI 13=C0S< X+PHl )*SIN< X+PHI ) * ( RHOQI *EKS“ 1 + 

1 3.  «COSX-«-ALOG(  (RHOQ<-RCGSX  ) / (EKS+RCOSX-1.  ))  + 

2 (4.  ■i^COSX^-J^O-l  ) / (EKS*<EKS+RHOQ-COSX  ) )- 

3 2.  s-COSX  f^  ( 1 . +ECOSX-RCOSX  ) / < EKG*  ( 1 . --C0SX^£”»-2  ) ) ) 

ELSE 

FI  13=C0S( X+PHI )*SIN(X+PHI )#<RH0QI*EKS-1.  + 

1 3.  ^^-COSX-«-ALOG( \EKS-RC0SX  + 1.  ) / ( RHOQ-RCOGX  ) ) + 

2 <4.  -t^C0SX«-fr2-i  > / (EKS*<EKS+RHOQ--CQSX  ) )- 

3 2 *COSXh^  ( 1 . -TCOSX  ) / < EKS*  ( 1 . -ECOSX-RCOSX  ) ) ) 

END  IF 

RETURN 

END 

FUNCriON  F023<X) 

COMMON  /PARFUN/  PHI, RHOQ, RHOQI, TRHOQ,  TRHOQI, RHOSQ,  RHOSQI, 03 
COSX=COS( X ) 

RC0SX=RH0Q*C03X 
TC0SX=2.  *RCOGX 
EKS=SQRT ( RHOSG-TCOSX ) 

EC0SX=EK3*C0GX 
IPUCDSX.  GT.  . 5)  THEN 

F023=SIN(  X+PHI  RHOQI  #EKB-1.  + 

1 3.  *C03X*AL0G( (RH0Q+RC0SX)/(EK3+RC03X-1.  ))+ 

2 (4,  ^^C0SX  t^H^2-  1 ) / (EKSXEKSh-RHOQ-COSX  ) )-- 

3 2 «-COSX-)t  < 1 . -TCOSX  ) / ( EKS*  ( 1 . -ECOSX-RCOGX  ) ) ) 

ELSE  IF(CGSX.  GT.  ~.  5)  THEN 

F023=BIN(  X+PHI  ) -«-«  2 a- < RHOQ  I -«<-EKS- 1 . + 

1 3.  -«-CGSX-«-ALOG<  (RH0Q'^RC0SX)/(EK5+RC0SX-1.  ))  + 

2 (4  ^^C0SXi^«2--  l.  ) / (EK5-«-<EKG-i-RH0Q“C0SX  ) )- 
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3 2,  -«-COSX»^(  1.  +EC0SX-RC0SX)/<EKS<^(1.  --C0SX*«2))) 

ELSE 

F023»SIN(  X4PHI  )*«2s-(RH0GI*EKS“l.  + 

1 3.  •i^COSX'Ji-ALOGv' (EKS-RCOSX  + 1.  ) / < RHOQ-RCOSX  ) ) + 

2 <4.  •»-C0SX**2-l.  )/<EKS»<EKS4-RH0Q-CDSX)  )- 

3 2.  «-CDSX*(l.  -TC□SX)/<EKG^^•(  1.  -ECOSX-RCGSX  ) ) ) 

END  IF 

RETURN 

END 

FUNCTION  F002(X) 

COMMON  /F^ARFUN/  PHI,  RHOQ,  RHOQI,  TRHOQ,  TRHOQI,  RHOSQ.  RHOSQI,  03 
SINX=SIM(X) 

COGX-^COS  ( X ) 

EKS2=RHDSGI~TRH0QI-«-C0SX 
I F ( AB5 ( S I NX ) . GE.  1 . E-4 ) THEN 

F02=ATAN(SINX/ (RHOQ-COSX) )/SINX 
ELSE 

DEN=1. /(RHOG-COSX) 

F02=DEN-03*S  I NX  *-fr2*DEN  »^-fr3 
END  IF 

F002=.  5« ALOG ( EKS2 ) +C0SX*F02 

RETURN 

END 


FUNCTION  F102(X) 

COMMON  /PARFUN/  PHI , RHOQ,  RHOQI , TRHOQ,  TRF^OQI , RH05G,  RHOSQI , 03 
SINX==--^SIN(  X ) 

COSX=-COS(  X ) 

EKS2=^RH0SQI--TRH0QI#C0SX 
IF(ABS(SINX).  GE.  l.E~4>  THEN 

F02=ATAN(SINX/  (RHOQ--COSX  > ) /SI MX 
ELBE 

DEM=1.  / <RF^0Q•~C0SX  ) 

F02=DEN-03->^S  I NX*«2*DEN-«-^^3 
END  IF 

F 1 02=^C0S  ( X +PH I ) * (RHOQ  I +COSX  »-AL0G  ( EKS2  ) 

1 - ( 1 . -2.  *CQSX-fr-R-2)*F02) 

RETURN 

END 

FUNCTION  F012(X) 

COMMON  /PARFUN/  PHI, RHOQ, RHOQI, TRHOQ,  TRHOQI, RHOSQ. RHOSQI,  03 
SINX==SIN(X  ) 

COSX=COS(X) 

C0SX2==C03X«-«-2 
EKS2=RH0SQI-TRH0QI«-C0SX 
IF(ABS(3INX).  GE.  I.  E-4>  THEN 

F02=ATAN(SINX/ (RHOQ-COSX) )/SINX 
ELSE 

DEN=1, /(RHOQ-COSX) 

F02=DEN-03*S  I NX*«  2^«-DEN^<-«  3 
END  IF 

FO 1 2=S IN  ( X +PHI  ) »^  ( RHOQI +COSX»ALOG  ( EKS2 ) 
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1 -(1.  -2.  t^C0SX**2)*F02) 

RETURN 

END 

FUNCTION  F202<X) 

COMMON  /PARFUN/  PHI, RHOQ, RHOQI, TRHOQ, TRHOQI, RHOSQ, RHOSQI, 03 
SINX:^SIN(  X ) 

COSX^=COS<  X ) 

C0SX2=C0SXi^i^2 
EKS2=Rh'0SQ  I -TRHOQI  *COSX 
IF(ABS(SINX).  GE  1.  E--4)  THEN 

F02=ATAN(5INX/ (RHOQ-COSX) )/SINX 
ELSE 

DEN==1.  /(RHOQ-COBX) 
t~02=DEN'~03'J^SI  NX‘fr*2*DEN'i^’'ti‘3 
END  IF 

F202=C0S(  X+PHI)  K-*2*(RH0QH<-( . 5*RH0QI+2.  *C0SX  ) 

1 +(2.  -!^C0SX2-  5)^^AL0G(EK52)-C0SX^^<3.  -4,  -t^C0SX2 ) *F02  ) 

RETURN 

END 

FUNCTION  F112(X) 

COMMON  /PARFUN/  PHI, RHOQ, RHOQI, TRHOQ, TRHOQI, RHOBQ, RHOSQI, 03 
SINX=SIN(X) 

COSX=COS( X ) 

C03X2=C0SX■fi■^^■2 
EKS2^--=RH0SQI-TRH0QI*C03X 
IF(ABG<SINX ).  GE.  1.  E-4)  THEN 

F02=ATAN(SINX/ (RHOQ-COSX ) ) /SINX 
ELSE 

DEN=1.  / (RHOQ  -CGSX  ) 

F02=DEN-  03  K-S I NX  2*DEN **3 
END  IF 

F112=C0S(X+-PHI  )*SIN(X-('PHI  )*(RHOQI*(.  5*RH0QI+2  ^tC0SX.) 

1 +(2.  *C0SX2-.  5)  »<-AL0G<EKS2)-C0SX-ii-(3.  -4.  *C0SX2)*F02) 

RETURN 

END 

FUNCTION  F022(X) 

COMMON  /PARFUN/  PHI, RHOQ, RHOQI.  TRHOQ,  TRHOQI,  RHOSQ,  RHOSQI, 03 
BINX=SIN( X ) 

C0SX=C0S(X) 

C0SX2=C0SX-»-«-2 

EKS2=RH0SQ  I -TRHOQ  I -k-COSX 

IF* ABS(SINX ) GE. l.E-4)  THEM 

F02=ATAN(SINX/ <RH0Q-C0SX ) ) /SINX 
ELSE 

DEN=1  /(RHOQ-COSX) 

F02=DEN- 03-«-S  I NX**2^H)EN>  fr3 
END  IF 

F022=SIN(  X+PHI  ) -K-*2-«-<  RHOQI  tt( . 5*RH0QI+2.  *C0SX  ) 

1 +(2.  -i^C0SX2-.  5)-«-AL0G<EKS2)-C0SX«(3.  -4.  *C0SX2)*F02) 

RETURN 

END 
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PROGRAM  SC ATP 

C 

C TO  COMPUTE  THE  SIZES  OF  ARRAYS  FOR  PROGRAM  SCAT  AND  RELATED  PROGRAMS 

C 

CHARACTER  QAL*4,  SUF*3,  JCAT«-2,  PR  I *2 
COMMON  /PARAM/  RAD> ALPHA, BETA, DT, DTMl . THDIR 

COMMON  /ICOM/  IX, IS, ISl. NTOT, Nl, N2, N3, NI, ISKIP, IPATCH, ICORR 
COMMON  /CONST/  C, CMl , P I , TOP  I 
DATA  C/3.  E8/ 

CM1  = 1.  /C 

PI  =4.  #ATAN(  1.  ) 

TOP  I •■=2.  *PT 
CONV,'=PI/lSO 

OPEM(  1,  FILE^-^  ^SCTQAL  STATUS^  'OLD  ' ) 

READ(  1,  #)  QAL..  ISKIP,  NT,  IPATCH,  PRI,  ICORR,  NTL,  NLPl 

PRINT  'INPUT  FILE;  SCIN'//QAL,  ' SKIP= ',  ISKIP,  ',  PATCH= ',  IPATCH. 

1 CORR=  ICORR 

□PEN  < 2,  F I LE-=^ ' SC  I N ' / /QAL,  STATUS^  ' OLD  ' ) 

READ (2,4  .)  DT,  RAD,  ALPHA,  BETA 
READ (2, 5)  NI 
READ  (2,  5)  N1,N2,  N3 
IF(N,t.  GT.  100)  STOP  1 
IF(N3.  EQ.  0)  N3=5 

C FIND  NUMBER  OF  PATCHES  AND  NUMBER  OF  NEIGHBORING  PATCHES 

CALL  RUNTl (LL12, LL20) 

NU-2.  «RAD-^CM1  «-DTMl  +6. 

I F ( N I . EG,  0 ) N 1=6.  »R AD/  ( C*DT ) 

IF(NI.LT. 0)  NU=1 

PRINT  'DT=',DT,  ',  RAD='.RAD,  '.  ALPHA= ',  ALPHA,  ' . BETA=',PE7A, 

1 NI=',NI,',  N1=',N1,',  M2=',N2,  ',  N3= ' , N3 

PRINT  -i^,  NU,  ' ROWS  NEEDED' 

C FIND  DIMENSION  OF  ARRAY  NEEDED  FOR  FAR-FIELD  COMPUTATIONS 

400  READ(2,  4,.  END=500)  THDIR,  PHDIR 
ITH=THDIR 
IPHI=PHDIR 
THDIR=THDIRi'rCOMV 
P HD I R = P HD I R ON V 
CALL  FARFL 

I F ( NTOT . GT , M I'M  ) NTM=NTGT 
GO  TO  400 
500  LL=I3+2 

C ADD  2 BECAUSE  COMPUTATIONS  ON  855  AND  205  SOMETIMES  DISAGREE 

LC=MAX0(NU,  NTM)4-2 

C PRELIMINARY  CALCULATIONS  ONLY 

IF(NI  LT.  0)  LC=1 

C FILE  WITH  DIRECTIVES  FOR  THE  FULL  SCREEN  EDITOR 

OPEN (5, FILE='SCTFSE', STATUS^ ' NEW ' ) 

LB=LL#*2/64+l 

LBl=LL/64+l 

C COMPUTE  NUMBER  OF  LARGE  PAGES 

NLP=  < LL* ( 2#LC+2*LL20+55 ) +20*LL 1 2+2* ( LB+LB 1 ) +65539 ) Z65536 
PRINT  *.  'LL=',LL,  ',  LC=',LC,  LL12= ',  LL12,  ',  LL20=',LL20 

PRINT  *,  'COMPUTED  NLP=',NLP 

IF(LL.  GT.  99999.  OR.  LC.  GT.  9999.  OR.  LL12.  GT.  999999.  OR  LL20  GT  99^) 
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1 STOP  10 

WRITE(5>  6)  LL,  LC,  LL12..  LL20 
C IF  NUMBER  OF  LARGE  PAGES  IS  GIVEN 

IF(NLP1.  NE.  O)  WLP=NLP1 
IF (NLP.  GT.  28.  OR.  NTL.  GT.  599940)  STOP 
C ASSIGN  TIME  LIMITS  (IF  NOT  GIVEN)  AND  PRIORITIES 

IF (NLP.  LE.  5)  THEN 

IF (NTL.  EQ.  O)  NTL  = 500 
IFdMTL.  LE.  900)  THEN 
JCAT= 'P2' 

ELSE  IF ( NTL.  LE.  3600)  THEN 
JC AT=  ^P3' 

ELSE 

JCAT='P5' 

END  IF 

ELSE  IF(NLP.  LE.  10)  THEN 
IF(NTL.  EQ.  O)  NTL=900 
IFdMTL.  LE.  3600)  THEN 

v.lCAT=  'P3' 

El.SE 

JC AT= 'P5 ' 

END  IF 

ELSE  IF (NLP.  LE.  20)  THEN 
IF (NTL.  EG.  0)  NTL=900 
IFdMTL.  LE.  900)  THEN 
JC AT= 'P4 ' 

ELSE 

JCAT=  'P5  •' 

END  IF 
ELSE 

IFdMTL.  EG.  0)  NTL=5000 
JCAT= 'P5 ' 

END  IF 

IF (PR  I.  NE.  ^XX  ' ) THEN 
IF (PRI. LT. JCAT)  THEN 

PRINT  "ASSIGNED  PRIORITY  TOO  LOW' 

PRI  NT  « . 'PRI  --=■  APR!.  ",  JCAT^^  ' , JCAT 
STOP 
ELSE 

JCAT^PRI 
END  IF 
END  IF 

0PEW(4, FILE= 'SCTRES ",  STATUS^ " NEW ' ) 

WRITE(4,  -K-)  'RESOURCE,  TL=",  NTL,  ",LP=",NLP,  " ..  JCAT=  " , JC  AT , 
1 ",  NT--- ",  NT,  ".  ' 

PRINT  «,  "PRIORITY  ",  JCAT,  ",  TIME  LIMIT  " , NTL , ",  ",NT.  ' 

STOP 

4 FORMAT (E 12  5) 

5 FORMAT (314) 

6 FORMATk  "REPLACE/LLG/ ",  15,  "./ALL  QUIET",/,  'REPLACE/LCG/ " 

1 "/ALL  QUIET ',/,  'REPLACE/LL12G/ ",  16,  "/ALL  QUIET",/, 

2 'REPL.ACE/LL20G/ ",  13,  "/ALL  QUIET  ",/,  "QUI T ' ) 

END 

r 


TAPE  " 
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SUBROUTINE  RUNT  I (LL12>  LL20) 


THIS  SUBROUTINE  DETERMINES  THE  NUMBER  OF  PATCHES  AND  THE  DISTRIBUTION 
NEIGHBORING  PATCHES  IF  ICORR  ~ 2 OR  3 

DIMENSION  NPHI(IOO) 

DIMENSION  XRK5000),  XR2<5000)>  XR3<5000),  THETA  < 5000) , PHI  I < 5000) 

COMMON  /PARAM/  RAD, ALPHA, BETA,  DTi  DTMl . THDI R 

COMMON  /ICOM/  IX,  IS,  ISl. NTOT, Nl, N2, N3,  NI,  ISKIP,  IPATCH,  ICORR 
COMMON  /CONST/  C, CMl , P I , TOP  I 

DT!-i=^P  I / ( N 1 - 1 ) 

BDTH2^SIN(DTHfr.  5) 

TH=PI 

IS=2 

NP)iJ  ( 1 

NPHI (Nl )=1 

N1H=(N1 +1 >/2 

DPH=-1. 

IF  ( IPATCH.  EG.  0)  THEN 
NPHX=-:.t 

90  NPHX=NPHX->i-2 

IFINPHX,  LT.  N3)  GO  TO  90 
DO  100  1=2, NIH 
TH=TH-DTH 
5TH=SIN(TH) 

NPI-N(N2-N3)-i-STH+N3+.  5 
IF(NPH.  GT.  NPHX)  NPHX  = NPHX-k-2 
NPHI ( I )=NPHX 
IS=IS4NPHX 

DPH=AMIN1 (DPH, STH*SI N ( P I /NPHX ) ) 

100  CONTINUE 

ELBE 

DO  105  1=2,  NIH 
TH=1H~DTH 
STH=3IN(TH) 

NP H=  ( N2  -N3  ) •i^-STH+  N3+ . 5 
NPHI  ( I .)=NPH 
IS=I'3+NPH 

DPH=AMIN1  (DPH,  STH-»^-SIN  ( P I /NPH ) ) 

105  CONTINUE 

END  IF 

I F ( MOD ( N 1 , 2 ) . EQ.  1 ) N 1 H=N 1 H- 1 
DO  1 1 0 J= I , N I - 1 
NPH=NPHI (NIH) 

NPHI ( J)=NPH 
I3=IS+NPH 
N1H=N1H-1 
110  CONTINUE 

DT1=1 . 6*RADKCM1*AMIN1 (SDTH2. DPH) 

IF(DT.  EQ.  0.  > THEN 
DT=DT1 
ELSE 

IF(DT1. LT. DT)  THEN 
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PRINT  3,  DT.DTl 
NI=NI«-DT/DT1 
DT=DT1 
END  IF 
END  IF 
DTM1=1. /DT 

PRINT  'NUMBER  OF  PATCHES  IS  IS 
C RETURN  IF  NO  NEIGHBORING-PATCH  CORRECTIONS  ARE  GOING  TO  BE  MADE 

IF(  ICORR.  LT.  2.  OR.  ISKIP.  EQ,  1 ) THEN 
LL12=1 
L!.20=l 
RETURN 
END  IF 

C CHECK  IF  ARRAYS  ARE  BIG  ENOUGH  TO  CALCULATE  THE  NUMBER  OF  NEIGHBORING 

C PATCHES 

IF( IS. GT.  5000)  STOP  5 

I3i==IS-NPH 

TH^^PI 

IN=1 

XRK  1 )~0. 

XR2(1 )=0. 

XR3<  1 )=-~RAD 
THETA ( 1 )=PI 
PHI  I ( 1 )=0. 

NP2=.l 

NP3=NPHI <2> 

INB2--=1 
INB3=2 
DPHI2=T0PI 
DPHI3=T0PI /NP3 
DO  200  I=2,N1-1 
TH=TH-DTH 
3TH=SIN(TH) 

CTH=COS(TH) 

RSTH==RAD-«-STH 
RCTH==RAD#CTH 
NPJ.=NP2 
MP2=NP3 
NP3=NPHI ( I+l ) 

INB1=TNB2 

INB2=INB3 

INB3=INB2+-NP2 

DPHI 1=DPHI2 

DPHI2=DPHI3 

DPH 1 3==T0'"  I /NP  3 

PHI=0 

DO  190  J=1,NP2 
IN=IN+l 
CPH=CGS<PHI ) 

SPH=SIN(PHI ) 

XRl  < IN)=RSTH-«CPH 
XR2( IN)=RSTHfrSPH 
XR3( IN)=RCTH 
THETA < IN)=TH 
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PHI  I ( IN)^PHI 
PHI=PHI+DPHI2 
IPO  CONTINUE 
200  CONTINUE 
XRl < IS)=0. 

XR2(  IS)=-"0. 

XR3(IS)«RAD 
THETA < IS) =0, 

PHII < IS)=0. 

DLl-~-DTH^fRAD 

MMAX==0 

IN=-^1 

NPNT^--0 

NP=0 

DO  240  1 = 1,  IG 
X=XR1 ( I ) 

Y=XR2( I ) 

Z=XR3( I ) 

TH=THETA( I ) 

PHI=PHI  Id) 

NM=  1 

:tF  (NP.EQ.  O)  THEN 
NPH=NPHI < IN) 

DPH=TOPI/NPM 
DP HP R = DP N*l.  1 
DL2=DPH*BI  N ( TH  > ->^R  AD 
RHAX=BQRTf2.  ) frAMAXl (DLl,  DL2 ) 

IFd,  EQ  1.  OR  I EQ.  IS)  RMAX=1.  5*DTH^^RAD 
NP=NPH 
IN^=IN+1 
END  IF 
NP=NP-1 
DO  220  J=1,IS 

DR=SQRT<  ( X-XRl  < J)  ) *#2-K  Y-XR2 ( J ) ) «-«-2+ ( 2-XR3 ( J ) )*«2) 

IF<  J.  EQ.  I.  OR.  DR.  GT.  RMAX)  GO  TO  220 
NN=NM+.l 
NPNT=NPM1 +1 
220  CONTINUE 

IF(NM.  GT.  NMAX ) NMAX=NM 
240  CONTINUE 

PRINT  *,  'NUNBER  OF  COLUMNS  NEEDED  IS  ',NMAX-1. 

1 ' TOTAL  NUMBER  OF  NEIGHBORING  PATCHES  IS  ',NPNT 

LI..12=NPNT 
I..L20=NMAX 
RETURN 

3 FORMAK'  GIVEN  DT=  ',E9.  2,  ' GREATER  THAN  "MINIMUM"  DT=',E9  2) 


SUBROUTINE  FARFL. 

C 

COMMON  /PARAM/  RAD, ALPHA, BETA, DT, DTM)  . THDIR 

COMMON  /ICOM/  IX, IS, ISl. NTOT, Nl, N2, N3. NI, ISKIP, IPATCH, ICORR 
COMMON  /CONST/  C, CM 1 , P I , TOP  I 
READ  < 2,1)  T1,T2 


SO 


C COMPUTE  FIRST  TIME.  OF  ARRIVAL  IF  T1  = 0 

IF(T1.  EQ  O.  ) T1=RAD*CM1*(  1.  -2.  K-SIN(.  S«^THDIR)  ) ~DT 
C COMPUTE  LAST  TIME  IF  T2  = 0 

IF<T2.  EQ.  0.  ) T2=(NI-5>*DT-RAD^<CM1 
READ (2, 2)  N 

MMAX=(TH-RAD*CM1  )/DT+4. 

NMIN=(T1--RAD*CM1  >/DT-3. 

NT0T=NMAX-NMIN<-1 
PRINT  'NTOT=MMTOT 
RETURN 

1 FORMAT (£12.  5) 

2 FORMAT (214) 

END 
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F'ROGRAM  PERF" 

C COMPUTES  MIE  SCArTERING  FOR  A PERFECT  CONDUCTOR 

PARAMETER  <Ny~a  192,  NTM=3000,  NT4==1200) 

CHARACTER  XLAB«-:30,  YLAB^<■30.  PTLABUtSO,  PTLAB2-K-80.  SBLAB*0O 
REAL  PT I( NTM ) , T I < NTM ) , AMP  1 < NTM ) , AMP2 ( NTM ) . WW ( NW ) , 

1  SI  < NW ) , S2  < NW ) , S3  ( NW ) , S4  ( NUI ) 

DIMENSION  PI (5, NT4 ) , P(NT4), TAU(D, NT4 ) , X(N14) 

COMPLEX  FT 1 < MW ) , FT2  < NW ) , FAC , EYE 
CGMMON/PERAR/ ANGL  E ( 5 ) , AG ,1  R ( 5,  NW ; , AG  1 U 5,  NW  ) , 

1 AG2R<5, NW); AG21 (3, MW) . 

2 BES  ( NT  4 ) ; EWEU  ( N 1 4 ) ..  R A ( N T 4 ) , £ I A ( NT4  ) , RB  ( NT4  ) , E I B ( NT4  ) .. 

3 P r [ , n ; AMP  1 , AMP2,  WW;  FT  1 ..  FT2;  S 1 , S2;  B3>  S4 , P I , P , TAU.  X 
COMMON/UAR/  IMANGL,  RHO,  EM 

DATA  XLAB/  'C-kTIME;  (METERS)  V>  YLAB/ ' IMTENSl  TY  V 
DATA  PTLABl/ 'SCATTERED  INTENSITY,  HORIZONTAL  POLARIZATION'/ 
DATA  PTLAB2/ 'SCATTERED  INTENSITY,  VERTICAL  POLARIZATION  '/ 
DATA  C/3.  Ee/,EYE/(0.  , 1.  )/ 

CMl^l. /C 
READ  H-,  I PLOT 

OPEN  ( 7,  FILE^  'PERFSV  '.  FORM=-  'UNFORMATTED  ' ) 

R.EAD(^«;  1 ) CTMAX,  CTMIN,  Rl,  ALPHA,  BETA 
READ ( *, 2 ) NT, NSTAR, NANGL 
C,DT=(CTMAX-CTMIN)  / (NT-1  ) 

IF (NANGL. GT.  5)  STOP  1 
1F(MSTAR.  GT.  MW)  STOP  8 
IF (NT. GT. NTM)  STOP  3 
READ  *,  (ANGLE( I ) , 1^1, NANGL) 

READ(-a  , 1)  DW,  WMAX 
NRE(.2=1 
CTIM=CTMIN 
DO  90  1 = 1,  NT 
PTI ( I )=CTIM 
TI ( I )=PTI ( I ) kCMI 
CTIM=CTIM+CDT 
90  CONTINUE 

CALL  ONEGA ( WW, WMAX , NGTAR , DW ) 

B2A2C2=  ( DE;:TA  k «-2-ALPHA*«2)  -KC*  fr2 
BAC  = - ( BETA-ALPHA  ) ^^C-«-2. 

ABC  = - ( AL  PHA+BETA  ) -«  C 
BC2=BETA-s^C-Jt*2 
ABC2=ALPHA«BC2 
MMAX=WMAX*Rl/C+7. 

IF(MMAX.  GT.  NT4-2)  THEN 
PRINT  'NMAX  ',MMAX 
STOP  20 
END  IF 
I x -=o 

TL=SECOND( ) 

DO  100  IW=2, NSTAR 
RH0=R1  )^WW(  IW)  / C 
EM=RH0+7. 

M=EM+.  5 

IF (M.  GT.  NT4-2)  STOP  21 
IX=IX+1 


tl2 


CALL  BESRJ(RHO,  . 5, M+2,  1, BES, NCALC) 

IF(MCALC.  NE.  Mh-2)  THEN 

PRINT  'ERROR  IN  COMPUTATION  OF  BESSEL  FUNCTIONS  ' 
STOP  10 
END  IF 

CALL  BESRY(RHO,  . 5,  l'H-2;  1,  ENEU.  NCALC) 

IF(NCALC.  NE.  h^2)  THEN 

PRINT  K-,  'ERROR  IN  COMPUTATION  OF  NEUMANN  FUNCTIONS  ' 
STOP  1 1 
END  IF 

C COMPUTATION  OF  SCATTERING  COEFFICIENTS 

CALL  SCACO 

C COMPUTATION  OF  MONOCHROMATIC  FIELDS 

CALL  FLDdW.  M) 

100  CONTINUE 

TL  = SECOND < )-TL 

PRINT  'TIME  TO  COMPUTE  MONOCHROMAT I C FIELDS  IS  ',TL,  ' 
DO  210  IAMG=1,NANGL 
FTl  ( 1 )=--0. 

FT2(  1 )=-0, 

DO  200  I--  2..  NSTAR 
W=WN ( I ) 

WRC=W-frRi/C 

FAC=CMPLX(B2A2C2..  EAC-i^W ) /CMPLX  ( ABC2- W t^^-2,  ABC*W;*^^2 
FAC=FAC-«-CMPLX  (COS(WRC  ) . SIN(WRC)  ) 

FAC  = BC2-»-FAC/Ui 

FTl  ( I )=EYE*CMPLX  ( AG1R<  lAlMG,  I ) , AGl  I < lANG,  I ) )rFAC 
FT2(  I )=^EVE^^-CMPLX  < AG2R  ( lANG,  I ).  AG2I  < lANG,  I ) )*FAC 
200  CONTINUE 

TL^SECGND( ) 

CALL  VINUFKTI,  AMP  1 , NT,  NSTAR,  WW,  1,  FTl,  SI,  S2,  S3,  S4,  -1  ) 
TL=SECOND( )-TL 

PRINT  K-,  'TIME  TO  COMPUTE  A FOURIER  TRANSFORM  IS  ',TL,  ' 
WRITE (7)  AMPl 

C COMPUTE  OTHER  POLARIZATION  UNLESS  THETA  =0,  ISO 

IF(AMGLE(IAMG).  NF.  0.  . AND.  ANGLEUANG).  ME.  180.  ) THEN 
TL=-S£COMf>(  ) 

CALL  VINUFTCTI,  AMP2,  NT,  NSTAR,  WW,  i,  FT2.-  SI,  52,  S3,  S4,  -1 
TL=-S£COWD(  )-TL 

PRINT  TIME  10  COMPUTE  A FOURIER  TRANSFORM  IS  ' , TL 
WRITE (7)  AMP 2 
END  IF 
210  CONTINUE 
REWIND  7 

DO  300  IANG-1,NAMGL 

WRITE(SBLAB,  3)  ALPHA,  BETA,  Rl>  ANGLE! lANG) 

IF(ANGLE( lANG).  NE.  0.  . AND.  ANGLE! lANG).  NE  ISO.  ) THEN 
READ! 7)  AMP 2 
READ! 7)  AMPl 
ELSE 

READ! 7)  AMPl 
END  IF 

DO  220  I^-1,NT 

AMPl  ( I )=AMP1  ! I )^<•■fr2 


, NCALC 
, NCALC 


SECONDS ' 


SECONDS  •' 

) 

. 'SECONDS  ' 
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220  CONTINUE 

IF<  IPLOT.  EQ.  1.  OR.  IPLOT.  EQ.  2)  THEN 

CALL  DRAW4< 1, 2, 16, 12, 44, 72, XLAB, YLAB, PTLABl, GBLAB ) 
CALL  DRAW4A(  1,  1,  NT,  ' ',  0,  0,  PTI,  AMPl,  0.  . 0.  ) 

CALL  DRAW4B<  1, 0,  0,  0,  NT,  PTI,  AliPl,  2.  , 2.  ) 

END  IF 

IF  (IPLOT.  EG.  0.  OR.  IPLOT.  EQ.  2)  THEN 

CALL  DRAW3( 1, 2,  16,  12, 44, 72, XLAB, YLAB, PTLABl,  SBLAB) 
CALL  DRAW3A( 1,  1, NT,  ' ^ 0, 0, PTI , AMP  1 . 0.  , 0.  ) 

CALL  DRAW3B(1,  0,  0,  0,  NT,  PTI,  AMPl,  2.  , 2.  ) 

END  IF 

IF < ANGLE (lANG) . NE.  0.  . AND,  ANGLE(IANG).  NE.  ISO.  ) THEN 
DO  260  1=1, NT 

AMP  2 ( I ) =AMP2  ( I ) •fi-s-2 
260  CONTINUE 

IF  ( IPLOT.  EG.  1.  OR.  IPLOT  EG.  2)  THEN 

CALL  DRAW4( 1, 2, 16, 12, 44, 72, XLAB, YLAB, PTLAB2, SBLAB) 
CALL  DRAW4A(1,  1,  NT,  ' ',  0,  0,  PTI,  AMP2,  0.  , 0.  ) 

CALL  DRAN4D<  1, 0,  0,  0.  NT,  PTI,  AMP2,  2 ,2.  ) 

END  IF 

IF(  IPLOT.  EQ.  0.  OR  IPLOT.  EQ.  2)  THEN 

CALL  DRAW3( 1, 2, 16, 12, 44, 72, XLAB, YLAB, P7LAB2, SBLAB) 
CALL  DRAW3A(1,  1,NT,  ' ' , 0,  O,  PT I , AMP2,  0.  ,0.  ) 

CALL  DRAW3B(  1, 0,  0,  0,  NT,  PTI,  AMP2,  2.  , 2.  ) 

END  IF 
END  IF 
300  CONTINUE 
CLOSE (7) 

CALL  G5PURG£<  O..FN= O^ERFSO',  ' STATUS=  ' , I STAT  ) 

IFdSTAT.  GT.  0.  AND.  ISTAT  NE.  1402)  STOP  4 
CALL  Q5DEFIME(  'LFN=  'PERFSU  ') 

STOP 

1 FORMAT (E 12.  5) 

2 FORMAT (14) 

3 FORMAT  ( 'ALPHA  = ',F0.  2..  BETA  = ',F  5.  2, 

1 RADIUS  U F4.  2,  ',  ANGLE  ',r4.0) 

END 

SUBROUTINE  SCACO 

PARAMETER  i NW  1 '"'2,  T!T  .QO-.-w,  m.  4=  1200) 

p ' I ‘ r.;!  ; . , T 1 ( NT  M > , AMP  1 ( I'n  ^ ^ AM: '2  ( im  i rl ; , NW  ( NO  ^ 
j.  b 1 ( NW  ) , S2  ( ) • ••"'3  ^ N', .: ) , S4  I'lN 
n T M'"t'.!C .[  OtM  F * ' NT4  ^ , P • i -i  i 4 ) , I AU  ( 5,  NT4  ) , X ( NT  4 ) 

COMPLEX  F 'l  1 (NN; , FT2(NW) , FAC,  EYE 
COMMON/PERAR /ANGLE ( 5 ) , AGIR ( 5, NW ) , AGl I ( 5, NW ) , 

1 AG2R ( 5, NW) , AG21 ( 5, NW) , 

2 BES(NT4), ENEU(NT4), RA(NT4), EIA(NT4), RD(NT4). EIB(NT4), 

3 PTI,  TI,  AMPl,  AMP2,  WW,  FTl,  FT2, SI , S2, S3, S4, P I , P. TAU,  X 
COMMON/ VAR/  NANGL, RHO, EM 

IERR=0 

M=EM 

C CALCULATION  FOR  PERFECTLY  CONDUCTING  PPUP-nr  r~ 

DO  82  I = 1 , M 
AY:=  I 


ATER=-(AY+.t.  )*e£S(  I + l )-RHO*BES(  1+2) 

BTER=(AY+1.  )*ENEU< I + l )-RHO*ENEU< 1+2) 

A2=ATER**2 

FRAC=1. /(A2+BTER**2) 

RB< I )=-A2*FRAC 

EIB ( I )=ATER*BTER*FRAC 

ATER-BES( I+l ) 

BTER=ENEU< I+l ) 

A2=ATER*#2 

FRAC=1.  /(A2+BTER**2) 

RA(I )=-A2*FRAC 
EIA( I )=ATER>BTER*FRAC 
82  CONTINUE 
RETURN 
END 
C 

SUBROUTINE  FLDdW,  LMAX) 

PARAMETER (NW=8192,  NTM=3000,  NT4=1200) 

REAL  PTI (NTM), TI (NTM), AMPl (NTM), AMP2(NTM),  WW(NW), 

1  SI (NW), S2<NW),  S3(NW),  S4(NW) 

DIMENSION  PI (5; NT4), P<NT4), TAU(5, NT4),  X<NT4) 

COMPLEX  FTl (NW), FT2(NW),  FAC,  EYE 
COMMON/PERAR/ANGLE(5),  AG1R<5,  NW),  AGII (5,  NW), 

1 AG2R(5,  NW),  AG2I  (5,  NW), 

2 BES(NT4), ENEU(NT4),  RA(NT4) , EIA(NT4) , RB(NT4),  EIB(NT4), 

3 PTI,  TI,  AMPl,  AMP2,  WW,  FTl,  FT2,  SI, S2, S3, S4, PI, P, TAU,  X 
COMMON/VAR/  NANGL,  RHO,  EM 

DO  3 J=l, NANGL 

THETA=ANGLE( J)*  1.  745329251994E-2 

P(l)=COS  (THETA) 

PI  ( J,  1 ) = 1. 

TAU( J, 1 )=P( 1 ) 

P(2)  = l.  5*P( 1 )**2-.  5 
PI  < J,  2)=3.  >P<  1 ) 

TAU(  J,  2)=4.  *P(2)-1. 

DO  2 L=3, LMAX 
EL=L 

P(L>  = (2.  -«-EL-l.  )/EL-«-P(  1 )*P(L-1  )-(EL-l.  ) /EL*P  ( L-2  ) 

PI ( J, L)=P( 1 )*PI ( J,  L-1 )+EL*P(L-l ) 

2 TAU( J, L)=EL*(EL+1.  ) *P < L ) -P < 1 ) *P I < J, L ) 

3 CONTINUE 

DO  4 L=l, LMAX 
EL=L 

4 X (L)  = (2.  *EL+1.  )/(EL*(EL+l.  )) 

M=EM+. 0001 

DO  6 J=  1 , NANGL 
AGIR  ( J,  IW)=0. 

AGII  ( J,  IW)=0. 

AG2R(  J,  IW)=0. 

AG2I  ( J,  IW)=0. 

DO  6 L=l, M 

AG1R<  J,  IW)=AG1R(  J,  IW)+X(L)*(RA<L)«-PI  < J,  L)+RB(L)  JtTAU<  J,  L)  ) 
AGIKJ,  IW)=AG1I(J,  IW)+X(L)*(EIA(L)*PI  ( J,  L)+EIB(L)«TAU(  J,  L)  ) 
AG2R(J,  IW)=AG2R(J,  I W ) +X ( L ) * ( RB ( L ) *P I ( J,  L ) +RA ( L ) *TAU ( J,  L ) ) 
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'w.' 


AG2KJ,  IW)=AG2I(  J,  IW)+X  <L)-k-<EIB<L)*PI  < J,  L)+EIA<L)*TAU(  J,  L)  ) 
6 CONTINUE  t . 

' ■)?';  I RETURN  ■■;■ » ; >••<;•>■  '*7  i 

END  ■ "' 
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PROGRAM  PERFPLT (INPUT, OUTPUT) 

C PLOTS  PREVIOUSLY  CALCULATED  SCATTERED  FIELDS 

CHARACTER  XLAD*30,  YLAB«30,  PTLAB1*80,  PTLAB2*80,  SBLAB^^SO 
DIMENSION  PTI (3000),  AMPOl (3000), AMP02(3000), 

1 T1 (300),  T2(300), T3(300), T4(300),  AMPl (300), AMP2(300), 

2 AMP3 ( 300 ) , AMP4 ( 300 ) , ANGLE ( 5 ) 

LOGICAL  Two, NEWT 1,  NEWT2 
CHARACTER ^<-4  (5L  1 , QLP..  QL3,  (5L4,  SL  ( 4 ) 

CHARACTER^s-40  L(6) 

DATA  XLAD/ 'C«TIME  (MEIERS)  '/, YLAB/ ' INTENSITY  '/ 

DATA  P7LAB1/ 'SCATTERED  INTENSITY,  HORIZONTAL  POLARIZATION'/ 
DATA  PTLAD2/ 'SCATTERED  INTENSITY,  VERTICAL  POLARIZATION  '/ 
DATA  C/3.  E8/, NEWTl, NEWT2/2*.  FALSE.  / 

READ  ■«,  GLl,  QL2,  QL3,  QL4,  CTl,  CT2,  IPLOT,  INVER.  (L(  I ),  1 = 1,  6) 

PRINT  •«,  'COMPARISON  OP  PLOTS  PRODUCED  BY  PERF  WITH  THOSE', 
i ' PRODUCED  BY  SCAT,  IF  ANY' 

PRINT  'QUALIFIERS  ' , QL 1 , QL2,  QL3,  QL4 
DO  100  1 = 1,6 
PRINT  4 
PRINT  *,  L(I) 

100  CONTINUE 
PRINT  4 

IF  (CTl.  NE.  0.  AND.  CT2.  NE.  O.  . AND.  CTl.  GT.  CT2)  THEN 

PRINT  'INITIAL  TIME  ',CT1.  ' GREATER  THAN  FINAL  TIME  ',CT2 
STOP  10 
END  IF 

IF  (CTl.  NE.  0.  ) NEWTl-=,  TRUE. 

IF(CT2.  NE.  0.  ) NEWT2=.  TRUE. 

CALL  Q5ATTACH( 'LFN=', 'PERFSV') 

0PEN(7, FIlE= 'PERFSV', FORM= 'UNFORMATTED ' , STATUG='OLD'  ) 

READ  1,  CTMAX, CTMIN. Rl, ALPHA,  BETA 
READ  2,  NT,  NSTAR..  NANGL 
C D T= ( C TM A X -C  TM I W ) / ( NT - 1 ) 

■ IF  (NANGL.  GT'.  5)  STOP  1 
IF (NT. GT.  3000)  STOP  3 
READ  ( ANGLE ( I ) , 1=1, NANGL ) 

READ  1,  DW, WMAX 

NREC=1 

NQ=0 

iF(QLl.  NE,  'MO'  ) THEN 
NQ=NQ4-1 
SL(NQ)=QL1 

CALL  Q5ATTACH(  •'LFN= ',  'SCFL'//QLD 
C)PEN(e,  FILE= 'SCFL '//QLl,  STATUS=  ' OLD  ' , 

1 ACCESS= 'DIRECT', FORM= 'UNFORMATTED ' , RECL=2400) 

END  IF 

IF( QL2. NE.  'NO ' ) THEN 
NQ=N(3M 
SL ( NQ ) =QL2 

CALL  Q5ATTACH(  'LFN= ',  'SCFL'//QL2) 

OPEN (9, FILE= 'SCFL '//QL2. STATUS= ' OLD ' , 

1 ACCESS= 'DIRECT ', FORM= ' UNFORMATTED ' , RECL=2400) 

END  IF 

IF(QL3.  ME.  'MO  ' ) THEN 
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NQ“NG+ 1 
SL ( NQ ) =QL3 

CALL  Q5ATTACH<  'i  'SCrL'//QL3) 

OPEN<  10,  FILE='SCFL  V/QL3,  GTATUG='□LD^ 

1 ACCESS^ 'DIRECT', F0RM= 'UNFORMATTED ' , RECL=2400) 

END  IF 

IF(QL4. NE.  'NO' > THEN 
NO-NGtl 
SL ( NQ ) =QL4 

CALL  Q5ATTACH(  'LFN= 'SCFL'//QL4) 

OPEN( 11, FILE= 'GCFL '//QL4,  STATUG= 'OLD 
1 ACCESS^-- 'DIRECT FGRM== 'UNFORMATTED  ' , RECL=2400) 

END  IF 

DD  500  IANG=1,NANGL 

WRrTE(SBLAB, 3)  ALPHA, BETA, Rl,  ANGLE! lANG),  ( GL < I ) , 1 = 1. NQ) 
IF(ANGLE(  lANG).  NE.  0.  . AND.  ANGLE(IANG).  NE.  180.  ) THEN 
TWO=.  TRUE. 

IF( INVER  EQ.  1 ) THEN 
READ (7, END=600 ) AMP02 
READ(7)  AMPOl 
ELSE 

READ ( 7, END=600 ) AMPOl 
READ(7)  AMP 02 
END  IF 
EL.SE 

TNO=  F-ALBE. 

READ ( 7 , END=600 ) AMPO 1 
ENDIF 

CTIM=CTMIN 
DC  210  1=1 , NT 

AMPO  id)  =AMP0 1 d ) -K-  5^2 
PTX ( I )=CTIM 
CTI  M=CTIM-iCDT 
210  CONTINUE 

I F ( NEWT  1 ) THEN 
DC)  222  101  = 1.  NT 

I F i P T 1 d 0 1 ) GE . C T 1 ) GO  TO  223 

222  CONTINUE 

PRINT  'INITIAL  TIME  OCTl,  ' TOO  LARGE' 

S T OP  1 0 
ELBE 
101  = 1 
END  IF 

223  IF(NEWT2)  THEN 

DO  224  I02~-NT..  1,  "1 

I F ( PT I < 1 02 ) . LE.  CT2 ) GO  TO  225 

224  CONTINUE 

PRINT  *,  'FINAL  TIME  ',CT2,  ' TOO  GMALL ' 

STOP  10 
ELBE 

I02=NT 
END  IF 

225  IF<I01.GT. 1)  THEN 

DO  226  13=101, 102 
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AMPOl  ( 13-101  -t-1  >=AMP01  ( 13) 

PTI ( I3-I01+1)=PTI ( 13) 

226  CONTINUE 

END  IF 

NTP=I02-I01+ 1 

IF<QL1.NE.  'MG')  THEN 
READ (8, END-600)  T1 
READ(8)  AMPl 
DO  240  NT 1=1. .300 

IF<  ANPl  (MTl  ).  LT.  0.  ) GO  TO  245 

240  CONTINUE 

245  NT1=NT1-1 

IF(NEWTl)  THEN 

DO  252  ,T1  = 1.NT1 

IF(T1 ( 1 1 ).  GE.  CTl ) GO  TO  253 

252  CONTINUE 
ELBE 

11  = 1 
END  IF' 

253  IF(NEWT2)  THEN 

DO  254  I2=1MT1,.  1 . -1 

IF(T1 ( 12) . LE.  CT2)  GO  TO  255 

254  CONTINUE 
ELSE 

I2=NT1 
END  IF 

255  IFdl.  NE  1 ) THEN 

DO  256  13=11, 12 

AMP  1 ( 1 3-1 1 +-1  ) =AMP  1(13) 

T1  < I3-I l-l )=T1 ( 13) 

256  CONTINUE 
END  IF 
NT1=12-I 1+1 

END  IF 

IF(QL2.  NE.  'NO ' ) THEN 
READ(9, END=600)  T2 
READ(9)  AMP 2 
DO  260  NT2=  1.-300 

IF(  AMP2(NT2).  LT.  0.  ) GO  TO  265 

260  CONTINUE 

265  NT2=NT2-1 

IF(NEWTl)  THEN 

DO  272  Ii=l,NT2 

IF(T2d  1 ) GE  CTl  ) GO  TO  273 

272  CONTINUE 
ELBE 

11  = 1 
END  IF 

273  IF(MEWT2)  THEN 

DO  274  I2=NT2> 1,-1 

IF(T2(I2).  LE.  CT2)  GO  TO  275 

274  CONTINUE 
ELSE 

I2=NT2 
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E£ND  IF 

275  TF< II.  NE.  1 ) THEN 

DO  276  1 3---=  I I,  12 

AMP  2 ( 1 3- 1 :l.  4 1 .)  ::*-AMP2  0. 3 .> 

T2(  I3“1 1 + 1 ):^T2(  13) 

276  CONTINUE 
END  IF 
NT2=I2-I 1+1 

END  IF 

IF(QL3.  ME.  'NOM  THEN 
READ<  10.  END==600)  T3 
RE AD (10)  AMP 3 
DO  200  NT3:=-1,300 

IF( AMP3(NT3) . LT.  0.  ) GO  TO  285 

230  CONTINUE 

285  NT3=NT3~1 

IF'INEWTl)  THEN 

DO  292  I 1=^1,  NT 3 

IF(T3< II ),  GE.  CTl ) GO  TO  293 

292  CONTINUE 
ELSE 

11  = 1 
END  IF 

293  IF(NEWT2)  THEN 

DO  294  I2=NT3, 1,-1 

IF(T3( 12).  LE.  CT2)  GO  TO  295 

294  CONTINUE 
ELBE 

I2=NT3 
END  IF 

295  IF( I 1 . NE.  1 ) THEN 

DO  296  13=11, 12 

AMP3( I3-I 1+1 )=AMP3( 13) 

T3(  I3-1 1 + 1 )-=T3(  13) 

296  CONTINUE 

END  IF 
NT3=I2-I  1 1 

END  IF' 

IF(0L4.  NE  'NO' ) 1 HEM 
R E AD ( 1 1 , END=600 ) T 4 
READ<il)  AMP4 
DO  300  NT 4 =1,300 

IF(AMP4(NT4) . LT.  0.  ) GO  TO  305 

300  CONTINUE 

305  NT4=NT4~1 

IF(NEWTl)  THEN 

DO  312  I1=1,NT4 

IF(T4(1 1 ) . GE.  CTl ) GO  TO  313 

312  CONTINUE 
ELBE 

I 1 = 1 
END  IF 

313  IF(NEWT2)  THEN 

DO  314  I2=NT4, 1, ~1 
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314 


315 


316 


320 


326 


340 
34  5 


IF(T4( 12).  LE.  CT2)  GO  TO  315 
CONTINUE 
ELSE 

I2=NT4 
END  IF 

IFdl.NE.  1)  THEN 
DO  316  I3=-1 1 , 12 

AMP4( 13- 1 1 + 1 )=AMP4( 13) 

T4 (1 3-1 1 + 1 )=T4( 13) 

CONTINUE 
END  IF 
NT4--=I2-1 1 + 1 
END  IF 
NREO-NREC+2 
NMl^^NT  1/20+4 
NM2-NM1 +1 
NM3^NN2+1 
NM4=NM3+1 

IFdFLOT.  EQ.  1.  OR  I PLOT.  EQ.  2)  THEN 

CALL  DRAW4( 1,  1,  16,  12,  44, 80,  XLAB,  YLAD, PTLABl,  SBLAB ) 

CALL  DRAW4Ad,  1..  NTP,  ' ' , 0,  0,  PTI,  AMPOl.-  0.  , 0.  ) 

IF(QL1.NE.  'NO')  CALL  DRAW4 A < 1 , 1 , NTl , ' 0 ' , NMl , 1 , T 1 , AMP  1 , 0.  , O.  ) 

IF(QL2.  NE.  3M0'  ) CALL  DRAW4Ad,  1,  NT2,  ' 1 ',  NM2,  2,  T2,  AMP2,  O.  , 0.  ) 

IF(QL3.  NE.  'NO')  CALL  DR  AW4A  d , 1 , NT3,  ' 2 ' , NM3,  3,  T3,  AMP3..  0.  , 0.  ) 

IF(QL4.  NE.  'NO')  CALL  DR  AW4  A d , 1 , NT4 . ' 3 ' , NM4,  4,  T4,  AMP4,  0.  , 0.  ) 

CALL  DRAW4Bd,  O,  0,  0,  3000,  PII,  AMPOl,  2.  , 2.  ) 

END  IF 

IF(  IPLOT.  EQ.  0.  OR.  IPLOT.  EQ.  2)  THEN 

CALL  DRAW3( 1 . 1,  16,  12, 44, 80,  XLAB, YLAB, PTLABl, SBLAB) 

CALL  DRAW3Ad,  1,  NTP,  ' ',  0,  0,  PTI , AMPOl , 0.  , 0.  ) 

IF(QL1  ME.  'NO' ) CALL  DR  AW3A  d , 1 , NT  1 ' 0 ' , NMl  , 1 , T 1 , AMP  1 , 0.  , 0.  ) 

IF(QL2  NE.  'NO')  CALL  DR AW3A  d , 1 , NT2,  ' 1 ' , MM2,  2,  T2,  AMP2,  0.  , 0.  ) 

IF(QL3.  NE.  'NO')  CALL  DRAW3A  d , 1 , NT3,  ' 2 ' , NM3,  3,  T3,  AMP3,  0.  , 0.  ) 

IF(QL4.  NE.  'NO')  CALL  DRAW3A  d , 1 , NT4,  ' 3 ' , NM4,  4..  T4,  AMP4,  O.  , 0.  ) 

CALL  DRAW3B  d,  0,  0,  0,  3000,  PTI,  AMPOl,  2.  , 2.  ) 

END  IF 
IF < TWO)  THEM 
DO  320  I==1,MT 

AMP  02  < I )-AMP02(  I )k-k-2 
CONTINUE 

IFdOl.NE.  1)  THEN 
DO  326  I '3--^^  1 01..  102 

AMP02(  1 3-101  + 1 )=^AMP02(  13) 

CONTINUE 
END  IF 

IF(QL1.  NE  'NOd  THEN 
READ (8, END-600)  T1 
READ(8)  AMPl 
DO  340  NT  1 = 1, 300 

IF(  AMPKNTl  ) . LT.  0.  ) GO  TO  345 
CONTINUE 
|\IT1=NT1-1 
IFHMEWTl)  IHEN 
DO  352  11=1, NTl 
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352 

353 

354 

355 

356 

360 

365 

372 

373 

374 

375 

376 


IF(T1  < 1 1 ).  GE.  CTl  ) GO  TO  353 
CONTINUE 
ELSE 
11  = 1 
END  IF 

IF(NEWT2)  THEN 

DO  354  I2=NT1,  1,  -1 

IF(T1 ( 12).  LE.  CT2)  GO  TO  355 
CONTINUE 
ELSE 

I2=NT1 
END  IF 

IF< II,  NE.  1 ) THEN 
DO  356  13=11, 12 

AMPK  I3-I1  + 1 )=AMP1  ( 13) 

T1  ( I3-1 14-1  )=T1  ( 13) 

CONTINUE 
END  IF 
NTl  = I2-Il4-i 
END  IF 

IF(QL2.  NE.  'NOM  THEN 
READ(9, END=600)  T2 
READ<9)  ANP2 
DO  360  NT2= 1,300 

IF(AMP2(NT2) . LT.  0.  ) GO  TO  365 
CONTINUE 
NT2=NT2-1 
IF ( NEWT 1)  THEN 
DO  372  11=1, NT2 

IF(T2< II ).  GE.  CTl ) GO  TO  373 
CONTINUE 
ELSE 
11  = 1 
END  IF 

IF(NEWT2)  THEN 

DO  374  I2=NT2,  1,-1 

IF(T2< 12),  LE.  CT2)  GO  TO  375 
CONTINUE 
ELSE 

I2=NT2 
END  IF 

IF( II.  NE.  1 ) THEN 
DO  376  13=11, 12 

AMP2<  I3-I14-1  )=AMP2(  13) 
72(13-11+1 )=T2( 13) 

CONTINUE 
END  IF 
NT2=I2-I 1+1 
END  IF 

IF(QL3.  NE.  'NOM  THEN 
READ (10, END=600 ) T3 
READ(IO)  AMP3 
DO  380  NT3= 1,300 

IF(AMP3(NT3).  LT.  0.  ) GO  TO  385 
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380 

385 

392 

393 

394 

395 

396 

400 

405 

412 

413 

414 

415 

416 


CONTINUE 
NT3=NT3-1 
IF ( NEWT 1)  THEN 
DO  392  11=1, NT3 

IF(T3(  II  ).  GE.  CTl  ) GO  TO  393 
CONTINUE 
ELSE 
11  = 1 
END  IF 

IF(NEWT2)  THEN 

DO  394  I2=NT3,  1,  -1 

IF(T3( 12).  LE.  CT2)  GO  TO  395 
CONTINUE 
ELSE 

I2=NT3 
END  IF 

IFdl.ME.  1)  THEN 
DO  396  13=11, 12 

AMP3( I3-I1+1 )=AMP3< 13) 

T3(  I3--1 1 + 1 )=T3(  13) 

CONTINUE 
END  IF 
NT3=I2-I1+1 
END  IF 

IF(QL4.  NE.  'HO' ) THEN 
READ( 1 1,  END=600)  T4 
READ(ll)  AMP4 
DO  400  NT4= 1,300 

IF(AMP4(NT4).  LT.  0.  ) GO  TO  405 
CONTINUE 
NT4=NT4~1 
IF ( NEWT 1)  THEN 
DO  412  11=1, NT4 

IF(T4< II ).  GE.  CTl ) GO  TO  413 
CONTINUE 
ELSE 
11  = 1 
END  IF 

IF(NEWT2)  THEN 

DO  414  I2=NT4,  1,  -1 

IF(T4< 12).  LE.  CT2)  GO  TO  415 
CONTINUE 
ELSE 

I2=NT4 
END  IF 

IF( II.  NE.  1 ) THEN 
DO  416  13=11, 12 

AMP4< I3-I1+1 )=AMP4( 13) 

T4( I3-I1+1 >=T4( 13) 

CONTINUE 
END  IF 
NT4=I2-I 1+1 
END  IF 
NREC=NREC+2 
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NMl«NTl/20+-4 

NM2«NM1+1 

NM3aNM2+l 

NM4«NM3+1 

IF<  IPLOT.  EQ.  1.  OR.  IPLOT.  EQ.  2)  THEN 

CALL  DRAW4( 1,  1,  16,  12, 44, 80, XLAB,  YLAB,  PTLAB2,  SBLAB) 

CALL  DRAW4A(  1,  1,  NTP,  ^ 0,  0,  PTI,  AMP02.  0.  , 0.  ) 

IF(QL1. NE.  'NO' ) CALL  DRAW4A< 1,  1, NTl,  ' 0 ' , NMl , 1 , T1 , AMP  1 , 0.  , 0.  ) 
IF(QL2.  NE.  'NO' ) CALL  DRAW4A(  1,  1,  NT2,  ' 1 ' , NM2,  2,  T2,  AMP2,  0.  , 0.  ) 

IF<QL3.  NE  'NO')  CALL  DRAW4A  < 1 , 1 , NT3,  ' 2 ' , NM3,  3,  T3,  AMP3,  0.  , 0.  ) 

IF(QL4.  NE  'NO')  CALL  DRAW4A ( 1 , 1 , NT4,  ' 3 ' , NM4, 4, T4, AMP4,  0 , 0.  ) 
CALL  DRAW4B  (1 , 0,  0,  0,  3000,  PTI , AMP02,  2.  , 2.  ) 

END  IF 

IF  < IPLOT.  EQ.  0 OR.  IPLOT.  EQ.  2)  THEN 

CALL  DRAW3< 1, 1, 16, 12, 44, 80, XLAB, YLAB, PTLAB2, SBLAB) 

CALL  DRAW3A (1 , 1 , NTP,  ' ' , 0, 0, PTI . AMP02,  O.  , O.  ) 

IF(QL1. NE.  'NO  ' ) CALL  DRAW3A( 1,  1, NTl,  ' 0 ' . NMl , 1 . T1 , AMP  1 . 0.  , 0.  ) 
IF(QL2.  NE.  'NO')  CALL  DRAW3A ( 1 , 1 , NT2,  ' 1 ' , NM2,  2,  T2.  AMP2,  0.  , 0.  ) 

IF(QL3.  NE.  'NO')  CALL  DRAW3A  ( 1 , 1 , NT3.  ' 2 ' , NM3,  3,  T3,  AMP3,  0 , 0.  ) 

IF(QL4.  NE.  'NO  ' ) CALL  DRAW3A ( 1 , 1 , NT4,  ' 3 ' , NM4, 4, T4, AMP4, 0.  , 0 ) 
CALL  DRAW3B(  1, 0,  0,  0,  3000,  PTI,  AMP02,  2.  . 2.  ) 

END  IF 
END  IF 

500  CONTINUE 
600  STOP 
1 FORMAT (E12.  5) 

FORMAT ( 14) 

FORMAT  ( ' ALPHA=  ' , F5.  2,  ' , BETA=  ' , F5.  2, 

1 ',  RADIUS= ',  F4.  2,  ',  ANGLE=  ' , F4.  O: 

2 QUAL.  ; ',A4;',  ',A4;',  ',A4:',  ' , A4 ) 

4 FORMAT <//) 

END 
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PROGRAM  CURCHK 

PARAMETER  <LL=LLG) 

C LCG,  LL12G,  AND  LL20G  NOT  NEEDED 

DIMENSION  CJTH(900,  1 0 ) , C JPH < 900,  10),T<900) 

DIMENSION  NPCH( 10), CJA(LL) 

CHARACTER  QAL*4, XLAB*30, YLAB*30,  TITLE*80,  SUBT*80 
COMMON  /PASS/  Nl, N2, N3,  IS,  IPATCH, RAD, C, DTI 
DATA  XLAB/'C  * TIME  V, YLAB/ 'SURFACE  DENSITY'/, 

1 TITLE/ 'SURFACE  CURRENT  DENSITY'/ 

DATA  C/3.  E8/ 

READ  INFORMATION  FOR  CURRENT  CHECK 

READ  *,  QAL, IPLDT, IPRINT, NN 
IF<NN.  GT.  10)  STOP  1 
READ  ■»•,  (NPCHd  ) , 1 = 1,  NN) 

READ  *,  CTl, CT2, IPATCH 

READ  IN  VALUES  OF  CONSTANTS  FOR  THE  PULSE 

READ  4,  DT, RAD, ALPHA,  BETA 
READ  5,  NI 
READ  5,  N1,N2,  N3 

CALL  RUNT2 
IF<DT.  EQ.  0.  ) DT=DT1 
IF(NI.EQ.  0)  NI=6.  *RAD/ (C-t^-DT) 

CDT=C*DT 
NI  l=CTl/CDT-f-l, 

NI2=CT2/CDT+1. 

CALL  Q5ATTACH(  'LFN=',  'SCSV'//QAD 

0PEN(2, FILE= 'SCSV '//QAL,  STATUS= 'OLD ',  ACCESS=  'DIRECT  ', 

1 FDRM= 'UNFORMATTED RECL=ISi{-S) 

IF(NI2.  EQ.  1 ) NI2=NI 
IF(NI1.  GE.  N12.  OR.  NI2.  GT.  900)  THEN 

PRINT  «,  'NI1=',  Nil,  ',  NI2= ',  NI2,  CT1=',CT1,',  CT2=',CT2 

STOP  2 
E.ND  IF 

NI I=NI2-NI 1+1 
NK=NI 1/12+2 
NK1=NK+1 
TT=(NI 1-1 )»CDT 
DO  200  I=NI 1, NI2 

READ<2,  REC  = I-k2-1  ) CJA(1;IS) 

DO  1 10  N=l, NN 

C JTH ( I , N ) =C JA ( NPCH ( N ) ) 
no  CONTINUE 

READ<2,  REC=I*2)  CJAddS) 

DO  120  N= 1 , NN 

C JPH  d , N ) =C JA ( NPCH ( N ) ) 

120  CONTINUE 

T( I )=TT 
TT=TT+CDT 
200  CONTINUE 

IF( IPRINT.  NE.  0)  THEN 
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PRINT  ^ INPUT  FILE  SCINSQAL 

PRINT  RAD=',RAD/'/  AL.PHA= ALPHA,  ' • BETA=',BETA, 

1 ^ MI=MMI,  S N1=',N1,S  N^=^N2,  M3~',N3 

DO  220  N= i , NN 
NU=NPCH<N) 

DO  210  I=NI1, NI2 

PRINT  2,  NU,  I,  CJTH(  I,  N),  CJPH(  I,  N) 

210  CONTINUE 

PRINT  3 

220  CONTINUE 

END  IF 

IF< IPRINT,  NE.  1 ) THEN 
DO  250  N=l, NN 

WRITE<SUBT, 1 ) QAL, NPCH(N) 

IFdPLOT.  NE.  1 ) THEN 

CALL  DRAW3( 1, 7, 8,  15, 23,  33,  XLAB, YLAB, TITLE, SUBT) 

CALL  DRAW3A<  1,  1,  NI  I,  ' 1 ' , NK,  0,  T ( NI 1 ) , CJTH(NI  1,  N),  0.  , 0 ) 
CALL  DRAW3A(1,  1,NII,  ' 2 ' , NKl , 0,  T ( NI  1 ) , C JPH< NI  1 , N ) , 0 ,0.  ) 
CALL  DRAW3B(  1,  0,  0,  0,  NI,  CJA,  T,  2.  , 2.  ) 

END  IF 

IFdPLOT.  NE.  0)  THEN 

CALL  DRAW4( 1, 7, 8,  15, 23, 33,  XLAB, YLAB,  TITLE,  BUBT) 

CALL  DRAW4A(  1,  1,  Nil,  ' 1 ' , NK,  0,  T < NI  1 ) , CJTH(NH,  N),  0.  , 0.  ) 
CALL  DRAW4A(1,  1, Nil,  ' 2 ' , NKl , 0, T < NI 1 ) , CJPH(NI1,  N),  0.  , 0 ) 
CALL  DRAW4B(  1,  0,  0,  0,  NI,  CJA,  T,  2.  , 2.  > 

END  IF 
250  CONTINUE 

END  IF 
STOP 

1 FORMAT  ( 'QUALIFIER  ',A4,  ',  PATCH  NUMBER  ',14) 

2 FORMAT( IX, 218, 2(3X,  1PE16.  9) ) 

3 FORMAT (///) 

4 FORMAT (E12.  5) 

5 F0RMAT(3I4) 

END 

SUBROUTINE  RUNT2 

C 

C THIS  SUBROUTINE  DETERMINES  THE  NUMBER  OF  PATCHES  AND 

DT  FROM  THE  PATCH  DISTRIBUTION  ON  THE  SPHERE 

COMMON  /PASS/  NI, N2. N3,  IS,  IPATCH, RAD, C, DTI 

C 

PI  =4.  *ATAN(  1.  ) 

DTH=PI/ (Nl~l ) 

SDTH2=SIN(DTH*.  5) 

TH=PI 

N1PN(N1+1 )/2 
DPH=1. 

IS=2 

IF ( IPATCH.  EO.  0)  THEN 
NPHX=1 

90  NPHX=NPHXit2 

IF(NPHX.  LT.  N3)  GO  TO  90 
DO  100  1=2, NIH 
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TH=TH-DTH 

STH=SIN(TH) 

NPH=<N2-N3)  i^STH+N3+.  5 
IF(NPH. GT. NPHX)  NPHX=NPHX*2 
IB=IS+NPHX*2 

DPH=AMIN1 (DPH,  STH*BIN(PI/NPHX ) ) 
CONTINUE 

IF(M0D<N1, 2) . EQ. 1 ) IS=IB~NPHX 
ELSE 

DO  105  1=2, NIH 
TH=TH-DTH 
STH=SIN(TH) 

NPH=(N2-N3)-*«-STH+N3+.  5 
IS=IS+NPH»2 

DPH=AMINi  <DPH,  STH-«SIN < P I /NPH  ) ) 
CONTINUE 

IF(M0D(N1, 2).  EQ.  1 ) IS=IS-NPH 
END  IF 

DT1  = 1.  6»^RAD/C#AMIN1  (SDTH2,  DFl!) 

RETURN 

END 
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SUBROUTINE  VNUFKX,  Y,  N.  NU>  0M>  lOM,  FT,  DX,  CX,  SX,  JFLAG ) 

C THIS  SUBROUTINE  CALCULATES  THE  FOURIER  TRANSFORM  OF  THE 
C FUNCTION  GIVEN  BY  STRAIGHT  LINES  JOINING  THE  POINTS  GIVEN  BY  THE 
C ARRAYS  X,Y.  THE  NUMBER  OF  IMPUT  POINTS  IS  N,  THE  FREQUENCY  ARRAY  OM  IS 
C PROVIDED  AND  HAS  DIMENSION  NU.  DX,  CX,  SX  ARE  REAL  ARRAYS  OF  DIMENSION  N 
C THAT  ARE  USED  FOR  SCRATCH.  THE  FOURIER  TRANSFORM  IS  GIVEN  BY  THE 
C COMPLEX  ARRAY  FT. 

C IF  IOM=l,  INPUT  AND  OUTPUT  OM  ARE  CIRCULAR  FREQUENCIES 
C IF  I0M=2,  INPUT  ARRAY  IS  CIRCULAR  FREQUENCY,  OUTPUT  IS  FREQUENCY 
C IF  I0M=3,  INPUT  ARRAY  IS  FREQUENCY,  OUTPUT  IS  CIRCULAR  FREQUENCY 
C IF  I0M=4..  INPUT  AND  OUTPUT  OM  ARE  FREQUENCIES 
C JFLAG^+l  IF  THE  FOURIER  TRANSFORM  HAS  A FACTOR  EXP(-HWT) 

C JFLAG=-1  IF  THE  FOURIER  TRANSFORM  HAS  A FACTOR  EXP(-IWT) 

DIMENSION  X(N),  Y(N),  DX(N.>,  CX(N),  SX(N),  OM(NU) 

COMPLEX  FT(NU) 

TOP  I =8.  *ATAN<  1.  ) 

TPIN=1. /TOPI 
IF( lOM.  LE.  2)  GO  TO  25 
OM (1 ; N ) -OM ( 1 i N ) *TOP I 
25  NP--N-1 

DX(  .t.;  NP)=QSVDELT<X(li  N)i  DX(1;  NP)  ) 

NO=----l 

IF(OM( 1 ) . EQ.  0.  ) THEN 
NO-2 

S-QSSSUM  < Q8VAD JM  < Y < 1 ; N ) ; NP ) *DX ( 1 ; NP ) ) 

FT( 1 )-CMPLX<S, O.  ) 

END  IF 

DX( 1; NP )=Q8VDELT( Y< 1; N); NP)/DX( 1;  NP ) 

DO  70  I-NO, NU 
W=OM( I ) 

WIN-1.  /W 

CX(  1;  N)=VCOS(W*X<  liN);CX(liN)) 

3X(1;  N)=--VSIN(W*X(  1;  N)i  SXd;  N)  ) 

S2R==Q8SSUM(Q8VDELT(CX(  li  N)i  NP)*DX<  1;  NP  ) ) 

S2I -QSSSUM ( Q8VDELT ( SX (1 ; N ) ; NP ) *DX ( 1 ; NP ) ) 

S1R=Y( 1 )«CX( 1 )~Y(N)*CX(N) 

SI  I-Y(  1 )itSX(  i ) -Y(N)-»J-SX  (N) 

FT(  I )-CMPLX  ( ( -SI  I +S2R-«-WIN)*WIN,  < S 1 R + S2I  *WI  N ) -«-W  I N ) 

70  CONTINUE 

IF(  lOM.  EQ.  1.  OR.  lOM.  EQ.  3)  GO  TO  85 
0M<  1;  N)-=OM<  l5  N)  ftTPlN 
85  IF( JFLAG.  EQ.  1 ) RETURN 

FT ( 1 ; NU ) -VCDNJG  < FT  < 1 i NU  > i FT  < 1 ; NU ) ) 

RETURN 

END 

C 

SUBROUTINE  VINUFT( X, Y, N, NU, OM,  lOM, FT, SCRR, SCRI,  CX,  SX,  JFLAG) 

C THIS  SUBROUTINE  CALCULATES  THE  INVERSE  FOURIER  TRANSFORM  OF  THE  FUNCTION 
C GIVEN  BY  STRAIGHT  LINES  JOINING  THE  POINTS  GIVEN  BY  THE  COMPLEX  ARRAY  FT 
C AS  A FUNCTION  OF  THE  REAL  ARRAY  OM.  THE  NUMBER  OF  INPUT  POINTS  IS  NU 
C THE  INVERSE  F.  T.  IS  ASSUMED  TO  BE  A REAL  FUNCTION  AND  GOES  IN  THE 
C ARRAY  Y,  AND  THE  REAL  ARRAY  X IS  GIVEN  AND  HAS  DIMENSION  N. 

C SCRR,  SCRI  ARE  REAL  ARRAYS  OF  DIMENSION  NU  AND  IS  USED  FOR  SCRATCH 
C IF  IOM-1,  INPUT  AND  OUTPUT  OM  ARE  CIRCULAR  FREQUENCIES 
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C IF  I0M=^2,  INPUT  ARRAY  IS  CIRCULAR  FREQUENCY,  OUTPUT  IS  FREQUENCY 

C IF  I0M=3,  INPUT  ARRAY  IS  FREQUENCY,  OUTPUT  IS  CIRCULAR  FREQUENCY 

C IF  I0M=4,  INPUT  AND  OUTPUT  ON  ARE  FREQUENCIES 

C vJFLAG=h-1  if  THE  INVERSE  FOURIER  TRANSFORM  HAS  A FACTOR  EXP(+IWT> 

C JFLAG=^-1  IF  THE  INVERSE  FOURIER  TRANSFORM  HAS  A FACTOR  EXP(-IWT) 

DIMENSION  X<N), Y<N), OM(NU), SCRR(NU),  SCRI (NU),  CX ( NU ) , SX < NU) 
COMPLEX  FT(NU) 

TOPICS.  *ATAN< 1.  ) 

TPIN=1. /TOPI 
PIN==2.  «-TPIN 

IF(vlFLAG.  EQ.  1 ) GO  TO  15 

FT (1 ; NU ) =VCON JG ( FT (1 ; NU ) ; FT ( 1 j NU ) ) 

15  IF<JDM.  LE.  2)  GO  TO  25 
OM ( 1 ; NU ) =OM  < 1 ; NU ) *TOP I 
25  MP=NU-1 

SCRI  ( 1.;  NP)  = 1.  /QaVDELT(OM(  li  NU);  NP  ) 

SCRR ( 1 ; NP ) =Q8VDELT  < VREAL  < FT  < 1 , NU ) ; NU ) ; NP ) *SCR I ( 1 ; NP ) 

SCRI ( 1; NP)=Q8VDELT<VAIMAG(FT( 1; NU), NU); NP)*SCRI ( 1 , NP ) 

DO  70  1=1, N 
W=X ( I ) 

IF(ABS(W).  LT.  1.  E-19)  THEN 

Y (I ) =P I N*Q8SSUM  < Q8VAD JM ( VREAL ( FT  < 1 ; NU ) ; NU  > ; NP ) * 

1 Q8VDELT<0M( 1; NU); NP) ) 

ELSE 

WIN=1. /W 

CX ( 1; NU)=VCGS( W*OM( 1; NU) ; CX  < 1 ; NU) ) 

SX ( 1 ; NU ) =VSI N ( W*OM ( 1 ; NU ) ; SX ( 1 ; NU ) ) 

S2R=QaSSUM(Q8VDELT(CX ( 1; NU); NP)*SCRR(1; NP ) - 
1 Q8VDELT  ( SX  ( 1 ; NU ) ; NP  ) <^SCR  I ( 1 ; NP  ) ) 

S1I=REAL(FT< 1 ) )*SX( 1 )+AIMAG(FT(l ) >*CXn ) 

1 ~REAL(FT<NU)  )-;^SX  (NU)-AIMAG(FT(NU)  )*CX  (NU) 

Y(  I )=PIN*(-S1 1-«-S2R*WIN)*WIN 
END  IF 
70  CONTINUE 

IFdOM.  EQ,  1.  OR.  lOM.  EQ.  3)  GO  TO  85 
0M( 1; NU)=OM< 1; NU)*TPIN 
35  IF( JFLAG.  EQ.  1 ) RETURN 

FT ( 1 ; NU ) =VCON JG ( FT ( 1 ; NU ) ; FT ( 1 ; NU ) ) 

RETURN 

END 
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